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ABSTRACT 


Two transient models have been developed to study the catalytic 
ignition in a monolithic catalytic reactor. The special feature in 
these models is the inclusion of thermal and species structures in 
the porous catalytic layer. There are many time scales involved in 
the catalytic ignition problem, and these two models are developed 
with different time scales. In the full transient model, the 
equations are non-dimensionalized by the shortest time scale (mass 
diffusion across the catalytic layer). It is therefore accurate but 
is computationally costly. In the energy-integral model, only the 
slowest process (solid heat-up) is taken as nonsteady. It is 
approximate but computationally efficient. 

In the computations performed, the catalyst is platinum and the 
reactants are rich mixtures of hydrogen and oxygen. One-step global 
chemical reaction rates are used for both gas-phase homogeneous 
reaction and catalytic heterogeneous reaction. The computed results 
reveal the transient ignition processes in detail, including the 
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structure variation with time in the reactive catalytic layer. An 
ignition map using reactor length and catalyst loading is 
constructed. 

The comparison of computed results between the two transient 
models verifies the applicability of the energy-integral model when 
the time is greater than the second largest time scale of the 
system. It also suggests that a proper combined use of the two 
models can catch all the transient phenomena while minimizing the 
computat ional cost . 
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CHAPTER I 


INTRODUCTION 


1.1 A Description of Catalytic Ignition in a Monolithic Reactor 

The catalytic combustor is a device which uses a catalytic 
active metal on the reactor surface to initiate and stabilize the 
combustion process. The granular bed catalytic igniter has been in 
use for many years in small monopropellant rocket thrusters for the 


altitude 

control of 

satel 1 i tes 

[ 1 ] . The 

use of 

monolithic 

bed 

catalytic 

combustors 

has been 

considered 

in gas 

turbines 

[2]. 


Recently, there has been interest in adopting a monolithic reactor 
for the ignition of hydrogen/oxygen mixture in rockets [3], This 
monolithic reactor has the potential of higher performance and 
longer catalyst life. A monolithic catalytic reactor consists of a 
solid substrate block with parallel channels (e.g. honeycomb). On 
the surface of the channel, a layer of porous material is coated. 
Within the pores, catalysts are deposited. One popular example is 
the use of platinum catalyst on alumina. 

In this work, ignition transient in a monolithic reactor is 
simulated. The model includes several parts : a non-reactive, 
impermeable solid substrate; a gas phase in the channel capable of 
homogeneous reaction; and in between, a porous catalytic layer where 
gaseous reactants can diffuse into and catalyt ical ly react. The 
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catalytic combustion processes in a monolithic reactor can be 
described in two phases: 

(I) Initially, a premixed fuel/oxidizer mixture is introduced into 
the reactor, gas reactants diffuse into the porous catalytic 
layer and are absorbed, a heterogeneous chemical reaction 
occurs with the help of the catalyst, and heat is released and 
transferred to the gas phase and substrate. At this stage, the 
gas-phase reaction is unimportant; it is catalytic-reaction 
controlled. After a short period of time, a sharp increase of 
heterogeneous reaction rate can be observed which 

characterizes the catalytic ignition. 

(II) Following phase I, the gas in the channel is gradually heated 
up by the catalytic layer, and the gas-phase reaction becomes 
active. Eventually, a sharp increase in the gas reaction rate 
occurs which characterizes the gas-phase ignition, and it ends 
phase II. Once the gas is ignited, the reactor goes to steady 
state rapidly, and all of the reactants are consumed completely 
short distance after the gas-phase ignition point. 


1.2 Literature Review 

Mathematical models applicable to a catalytic combustor can be 
divided into two categories ; steady-state models and transient 
models. The most comprehensive steady-state model was one 
developed by Kelly, et al. [4]. It is one-dimensional and includes 
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homogeneous and heterogeneous reactions. In addition, the solid 
heat conduction, radiation and axial diffusion are considered. 
Their predictions show that the homogeneous reactions in catalytic 
combustors are " flame-like " and that beyond a certain mass flow 
rate they are extinguished. The model successfully predicts the 
effect of design parameters such as channel diameter and inlet gas 
temperature on operating characteristics such as blowout and 
breakthrough. For transient models, T’ ien [5] was the first one to 
deal with a monolithic reactor. He developed a quasi-steady, 
one-dimensional model in which the substrate energy equation is 
assumed to be transient while the gas-phase energy and species 
equations are the steady ones. This assumption is based on the fact 
that the longest characteristic time scale in the monolithic 
reactor, because of its nature large thermal inertia, is the 
substrate heat-up time. In addition, heat conduction along flow 
direction in the gas channel and the solid are neglected for high 
speed flow. This model correctly predicts transient behavior and is 
highly efficient in computation. Sinha, et al. [6] proposed a 
complete coupled two-dimensional transient equation for both the gas 
and solid phases. But the computation time was too long for 
parametric studies and the quasi-steady gas phase approximation was 
used for practical purpose. Prasad, et al. [7] developed a 
one-dimensional transient model incorporating axial conduction and 
heat losses but neglecting radiation. For small times, the model is 
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fully transient and solves time-dependent equations for solid 
temperature, gas temperature, and gas species concentration. 

However, because of the small time step needed to solve these 
equations, making the solution computationally expensive. a 
quasi-steady model is used as in [5], One of findings using this 
model is that the length required for achieving complete conversion 
in steady-state operation is shorter than that required in transient 
operation. Kesten, et al. [8] developed a two-dimensional full 
transient model which permits computation of concentration and 
temperature in the bulk gas phase and within porous catalyst 
particles. This model was applied in a reactor packed with porous 
catalyst particles; the reactants were rich-hydrogen/oxygen 

mixtures. 

Much of the modeling work on the monolithic catalytic combustor 
has assumed that heterogeneous reactions occurred on the surface of 
catalytic layer only. In reality. the catalyst particles are 
deposited in the pores and during reaction there can be a 
concentration and thermal structure within that layer. Problems 
involved concentration and temperature transients in a catalyst 
pellet [9,10] which, used extensively in chemical reactors and 
automobile converters, has been studied experimentally and 
theoretically in chemical engineering literature. These results 
have shown the existence of temperature and concentration gradients 
inside the pellet because of porous structure effect. They suggest 
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that an improved monolithic catalytic combustion model can be made 
by incorporating a porous catalytic layer with finite thickness. A 
more detailed literature review for catalytic combustion is offered 
by Pf efferle and Pfefferle [11]. 

1.3 Present Model 

Two transient models have been developed in this work to study 
the ignition process in a monolithic catalytic reactor. A rich 
hydrogen/oxygen mixture is used as the reactant, and the catalyst is 
platinum. One-step global reaction rates are used in both the 
homogeneous and the heterogeneous reactions. These models differ 
from the previous ones in that they include the effect of a finite 
porous catalytic layer, where reactants’ diffusion, consumption and 
heat conduction are important for the results of ignition phenomena. 

Two models are developed because of the consideration of 
computational time. In the first model, all transient processes are 
included. It is therefore a more complete model but is 
computationally costly. A second model is formulated with only the 
longest transient (solid heat-up) included in the unsteady term. It 
is more computational efficient and is a direct extension of the 
work in [5] . 

In chapter II, these models are formulated. Numerical methods 
are presented in Chapter III. Computed results using both models 
are given in Chapter IV, both for clarifying the transient ignition 
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process and for comparing the predictions from two models 



CHAPTER II 


MODEL DEVELOPMENT 


2.1 Introduction 

A monolithic reactor is composed of many uniform parallel 
channels; if the reactor is well insulated, then the study of the 
whole reactor can be reduced to the consideration of a single cell 
unit. Each cell consists of an open gas channel and an associated 
solid fraction (Fig. la). The solid is composed of two parts : a 

porous catalytic layer where catalysts are deposited in the pores 
(Fig. lb) and an impermeable solid substrate which provides monolith 
support. Commercially, the porous layer is often referred to as the 
washcoat because the catalysts are distributed inside the pores, 
reactants have to diffuse into the pore before being absorbed on the 
surface and reacted catalytically. The rate of in-depth reaction 
can therefore be regulated by the diffusion rate inside the pore. 
On the other hand, although the substrate does not directly 

participate in the reaction, during transient ignition, its large 
thermal inertia can have a great influence on the catalyst 
temperature and the catalytic reaction rate. The heat transfer 

through the solid and its heat-up process therefore have to be 
included in the ignition model. Figure lc specifies the symbol 
designations for the dimensions of gas channel, catalytic and 
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substrate layers. Although the model to be presented can be a 
variety of transient catalytic combustion processes, the particular 
calculations performed are for the start-up ignition transient of a 
rich-hydrogen/oxygen mixture. It is assumed in the theoretical 
calculation that the reactor is filled with hydrogen; then the 
deficient species oxygen is admitted at a rate which yields a 
specified mixture mass ratio. The oxygen convects downstream 
through the gas channel. In the meantime, it diffuses into the 

pores in the catalytic layer and becomes absorbed. Catalytic 

reaction occurs with a rate depending on the local oxygen 
concentration and temperature. The heat released from the reaction 
is conducted through the catalytic layer to both the substrate layer 
and the gas in the channel. As the gas temperature is raised 

sufficiently, the gas-phase reaction is initiated. This process 

quickly increases the gas temperature further to that of the 
adiabatic flame temperature and the ignition is completed. 

It is clear that transient phenomena in the monolithic 
catalytic combustor are complex due to many time scales involved. 
Following T’ ien [5], Table 1 defines these characteristic times and 
their estimated magnitudes. Several new time scales, which account 
for the in-depth heat/mass transfer processes in the catalytic 
layer, are added or modified. Added are the heat and mass transfer 
scales in the porous layer. The heterogeneous reaction, previously 
regarded as a surface reaction, is now treated as a volumetric 
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reaction. 

An one-step overall heterogeneous reaction rate in the 
catalytic layer is adopted from available experimental data [12-14]. 
The experimental results show that in a hydrogen-rich mixture, the 
reaction rate depends on the concentration of deficient oxygen only. 
Therefore, only the species for oxygen is required in the model. 
The derivation of the catalytic reaction rate expression can be 
found in Appendix A. 3. For the gas-phase reaction, an one-step 
overall rate is also assumed. The ignition delay data [15-17], 
matched separately with the expression from the spontaneous ignition 
theory [18], allows the rate constants to be determined (see 
Appendix A. 2). Although the elementary reaction kinetics for H 2 /O 2 
mixture are known, we have adopted the overall rate expression in 
order to be consistent with the heterogeneous reaction (for which 
elementary reactions are not fully understood). 

2.2 Governing Equations 

The complete set of equations needed to describe the transient 
behavior include that (1) the continuity, energy, and species 
equations in the gas phase; and (2) the energy, species equations in 
the solid. The following assumptions are made in the derivation of 
these equations : 

(a) The plug flow approximation is used in the gas channel and the 
pressure is assumed to be uniform. 
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(b) The heat/mass transfer coefficients between the gas and the 
solid are functions of distance and temperature, and derived 
from Nusselt number correlation in entrance flow in a tube 
[19]. A local stagnation flow estimation [4] is made for the 
first point in the entrance flow. 

(c) Heat conduction and mass diffusion in the axial direction are 
neglected because the Peclet number based on typical gas 
velocity O 10 M/S) is much greater than unity [20]. 

(d) Radiation heat transfer in both gas and solid phases are 
neglected. 

(e) An effective diffusion coefficient is used to calculate oxygen 
diffusion through the porous media. The effective diffusivity 
is assumed to be proportional to the reciprocal of local gas 
temperature only, i.e. P^> e = const. [21]. 

(f) The specific heat of mixture is constant. 

Gas-phase 

The general dimensional equations are listed as follows. 

equation of state P = P R (2.1) 

dp * 9 (p*u ) 

.. .. + = 0 ( 2 - 2 ) 

continuity ; + • 

at a x 
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momentum 


P = P (x=0 , t=0 ) 


(2.3) 


d Y 

* g * « 

species p ~ + p u 

3 t* 


d Y h S 

g d * 

; + — (Y - Y )= W 

3 X A g c,s g 


(2.4) 


. . 3 T „ ...ST* h* S* 


energy p C 2 +p u c § + — L_ ( T * _ T * )= q * w * 

p -- D * * g C,s n g 


P d t P d x A 


(2.5) 


where the reaction rate is given by (see Appendix A. 3 for results): 


-E 

W = K P a p b exp( 2-) 

g on * * 

9 R T 

g 


( 2 . 6 ) 


Defining nondimensional variables by 


x = 


t = 


q = — — 

a * * 

C T(o,o) 

p 


(2.7) 


u = 


U (0,0) 


Y = 


og 


y = 


g Y (o,o) 
g 


P = — 


p (0,0) 


T = 
g 


T(o,o) 


E = 
g 


R T (o,o) 
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The dimensionless equations are : 


dp d (pu) = Q 

continuity % 3 t d x 


( 2 . 8 ) 


species 


energy 


3 y 


a y 


oS — - + pu -r — — + J (y ~ y ) 
p °q at H a x d g c,s 

-E 

= B* (py ) b exp ( — =?— ) 
i <3 l „ 


a t 


a t 


o5 2 + pu + J (T - T ) 

P g at h 3x H g C,s 

-E 

= q B* (py ) exp ( T g - ) 

9 2 K g T g 


equation of state 


P 


= p T 


g 


(2.9) 


( 2 . 10 ) 


( 2 . 11 ) 


where the length scale in x-direction is chosen to be the reactor 
length L*. the time scale t* used for nondimensionalization is yet 

to be determined, and 


L / u (o.o) (2.12) 

V — 

The system of gas-phase equations (2.8) to (2.11) needs initial 
conditions. In this study, initially the reactor is filled with 
hydrogen; the start-up is from a step-wise oxygen injection. 
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Initial conditions 


<< 

\Q 

X 

II 

o 

II 

O 

II 

T (x, 

<3 

II 

o 

II 

4-J 

(2. 13) 

p (x, t=0) = 1 ; 

u (x, 

t=0) = 1 



Solid-phase 

Solid phase consists of two layers, a catalytic layer on top of a 
substrate layer. The dimensional governing equations are : 


species in catalytic layer 


€ r (- 


* a P 


* 

oc 


a t 


) = 


a r 


. . . 8 Y c . . . a 2 v 

- (r p D — £ ) + r p D - r V 


a r 


[2. 14) 


a x 


* 

where is the effective oxygen diffusivity and D* <x 


Tg 


In 


this work, we choose y-1 and it gives a constant value for the 
product p* D*. Thus Eq. (2.14) becomes 


‘-H? > r '^ »* h, } 


(2. 14a) 
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Boundary conditions 


» » 

at r = r 

d 


* 

at r = 


* 

r 

c 


h* p*(Y - Y ) 

d c g 



8 Y 

c 

* 

3 r 


3 Y 

c 

* 

3 r 


(2. 15a) 


(2. 15b) 


energy in catalytic layer 


^ C 


3T r « 3 T 3T\ * * 

*( | ) - A* { (r* 1 ) - ~^ 2 } * q W c 

c a t c l r a r a r a X -> 


( 2 . 16 ; 


Boundary conditions 


* * 

at r = r 

d 


a t 


* * ^ , c \ 

h (T - T ) = A ( — ) 


t c g 


a r 


(2. 17a) 


* » 

at r = r 

C 


a t a t 

A* ( J-) = A ( ?-) 

c a r a r 


(2. 17b) 


reaction rate 


t * 0.8 -E 

W* = 0*1.772 (p* T*) exp( — ^ ) 

C OC c R J 


(2. 18) 


where 0* is the catalyst loading. It describes the amount of 
catalysts deposited in the porous layer and the unit is [kg-pt/m 3 ] 
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where it is based on the substrate volume. Equation 
derived in Appendix A. 

energy in substrate layer 


* * 

P C ( 
s s 


a t 


d t 


* ( * * d T a T x 

- * | -^-r <r* _i ) * _f l 

S ^ r a r a r* d x* 2 / 


Boundary conditions 


at r = r 


A (- 

S 


d T 


d r 


-) = A (- 


a t 


a r 


* * 

at r = r 

s 


a t 


a r 


= 0 


Before nondimensionalizing the solid equations, the 
quantities are defined : 


* 

a = 

c 


* * 

P C 

c c 


a - 
s 


* * 

P C 

S s 


*2 


* 

T = 
t 


a 


€ h 


*2 


T = 


D (o,o) 

e 


Sh = 


* * 

D /h 

e c 


Be = 


A /h 

c c 


(2.18) is 


( 2 . 19 ) 


(2.20a) 


(2.20b) 


following 


( 2 . 21 ) 
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K 


* 

c 


M -,o. 8 

8* 1.772 p* Y (0,0) —f- exp(-E ) 

g M -I 


D = 

at 


* 

K 

c 







D 

at 


q 


c 


D = 

am 


K* h 


*2 

c 


p (o,o)D ( 0 , 0 ) 

e 


/ 

D = 

am 


D 

am 

Y (0,0) 
<? 


where a* and a* are the thermal diffusivity in the catalytic and 

C s 

substrate layer, respectively, x* is the thermal diffusion time in 
the catalytic layer (h‘ is the catalytic layer thickness), x* is the 
mass diffusion time in the catalytic layer, Sh and Be are the 
Sherwood and Biot number, respectively, K* is a group number from 
the catalytic reaction rate, and are the Damkohler 

numbers which have the following physical meaning : 


catalytic reaction rate 
0 mass diffusion rate in catalytic layer 


catalytic reaction rate 


D « 

at 


thermal diffusion rate in catalytic layer 


(2.22b) 
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Several time ratios are also defined : 


6 = 

c 



5 = 

s 


*2 * 

h / a 

c s 


T 


(2.23) 


5 = 

m 



Introducing the nondimensional parameters : 


r = 





(2.24) 


Y = 

c 



E = 

c 


E 


c 


R T (o,o) 


q = 

c 


* * 


C T (o,o) 

c g 


y = 

c 


Y 


Y ( 0 , 0 ) 
g 


T = 

c 


T 

C 

T (o,o) 
g 


A = 



Eq. (2.14a) for 02 becomes : 


5 

m 


T 


g 


d y 

c 

ITT 


r a r 


(r- 


5 y 


— ) - D 

r am 


0 . 8 

y_ exp 


E (1- 

C 


1 

T 


(2.25) 



18 


Boundary conditions 


at r = r 

d 


a y 

c 

d r 


Sh (y - y ) 

c g 


(2.26a) 


at r = r 

C 


a y 

c 

aT 


= o 


(2.26b) 


energy in catalytic layer 


5 

c 


a t 
ITT 


a 

r a r 


(r 


a t 
TT 


+ 


' 0. 8 
D y exp 

at c 



(2.27) 


Boundary conditions 


at r = r 

d 


a t 
TT 


Be (T - T ) 

c g 


(2.28a) 


at r = r 

C 


a T d T 

— = * —s — ~ 

3 r 3 r 


(2.28b) 


energy in substrate layer 


5 

S 


a t 


a t 


r 


a t 
TT 


) 


r a r 


(2.29) 
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Boundary conditions 


at r = r 

c 


d T d T 

dr dr 


(2.30a) 


at r = r 

S 


d T 

ITT 


0 . 


(2.30b) 


Initial condit ions 
catalytic layer 

T (x, r, t=0) = 1 

c 

substrate layer 

T (x, r, t=0) = 1 


y c (x, r, t=0) =0 (2.31) 


(2.32) 


In deducing the above system of equations, it is assumed that 

* * 

h c /L « 1. Hence the diffusion and conduction in the x-direction are 
negligible compared with those in the r-direction. It is assumed 
that we have a cold start-up initial condition, i.e. the solid 
temperature is uniform and equals the inlet gas temperature at the 
initial moment (t=0), and the catalytic layer has no absorbed 


oxygen. 
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2.3 Choice of Reference Time Scales 

In the dimensionless equations, a number of time ratios appears 
but the characteristic time scale has yet to be specified. The 
choice of the proper reference time scale is clearly dictated by the 
physical processes we would like to investigate. In Table 1, 
various time scales are estimated. One sees that mass diffusion 

time in the catalytic layer is the shortest while the solid heat-up 
time is the longest. If we choose mass diffusion as the time scale 
to be resolved, then we clearly have accounted for all the 
transient events in the combustor (except for acoustic wave 
phenomena). Such a model will be referred to as the " full 
transient model " . The full transient model is, in principle, the 
most accurate. But the calculation can become very time-consuming 
if it is carried out over a period of solid heat-up time. This is 
because the ratio of solid heat-up time to mass diffusion time is in 
general a very large number. So the full transient model is best 
reserved for instances where a fast transient events need to be 
reso Ived . 

On the other hand, if we choose the longest solid heat-up time 
as the reference time scale, then Table 1 shows that all other 
processes can be regarded as quasi-steady. Although fast transient 
phenomena can not be captured, the model can describe the global 
heat-up/ignition process and is fast in computation. Such a model 


will be referred to as the 


energy integral model 


In the 
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following section, both the full transient and the 
model will be derived or formulated. 


2.4 Full Transien t Model 

Substituting x by x* into Eq. (2.23), we get 

m & 


5 = 
g 


L / u (o,o) 


T 

m 


5 = 1 

m 


5 = 

c 


T 

m 


5 = 

s 


*2 * 

h / a 

c s 


T 

m 


and Eqs. (2.8, 2.9, 2.10, 2.11, 2.25, 2.27, 2.29) are 


Gas-phase 


continuity 


6 

g 


d P 
d t 


+ 


d (pu) 
d x 


species 


d y d y 

ps_ -^3. + pu lr -3 + J (y - y ) 


g 5 t 


v n / J 

OX D g C , s 


-E 


= B (py ) exp ( — =2- ) 
1 g r 


energy-integral 


(2.33) 


reduced to 


(2.34) 


(2.35) 
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energy 


p5 


d T 

q 

d t 


pu 


a t 

9 

a x 


+ j 

H 


(T - T ) 
g c, s 


-E 

= q B (py ) exp ( — =^- ) 
g 2 g i 


( 2 . 36 ) 


equation of state 


P = p T 

g 


( 2 . 37 ) 


Solid phase 


species in catalytic layer 


T 


g 


a y 

c 

a - T 


r a r 


(r- 


a y 

c 

~d~r 


' 0.8 
D y exp 

am c 




( 2 . 38 ) 


energy in catalytic layer 


5 

c 


a t 
ITT 


r a r 


(r 


a t 
TT 


+ d' 

a t 


0 . 8 

y exp 

c 



( 2 . 39 ) 


energy in substrate layer 


S 

S 


a t 


a t 


a 

r a r 


a t 


a r 


( 2 . 40 ) 
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2.5 Energy Integr al Model 

Consider a control volume in the solid phase which includes the 
catalytic and substrate layer (see Fig. lc). Three parts 
constitute the overall transient solid energy balance : 

(1) Time rate of change in internal energy in the control volume 
which includes the catalytic layer from r* to r* and the 

d c 

* * 

substrate layer from r to r . 

c s 


# 2 p C ti r 


J r d 


a t 

* 

d t 


* * 

dr dx 


(2.41) 


(2) Time rate of energy transfer into the control volume from the 
gas phase 


h 2 7T r dx (T* - T* ) (2.42) 

t d g c, s 


(3) Energy generated in the control volume per unit time 


r r 

s 

* * * * * 

W q 2 n r dr dx (2.43) 

w C 

r 

J d 


Taking energy balance, i.e. (1) = (2) + (3), the result is: 


* * a T 

m P C r dr = h r.(T 


a t 


t d g 


T ) + 

c, s 


W q r dr 

* c 

r 

d 


(2.44) 
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Eq. (2.44) can be used to define a solid heat up time. By equating 
the unsteady term with the convection term, we get 


* (i ) 
r 

HT 


pc r dr 

* n 


(2.45: 


Nu X (0,0) 
oo g 


If we equate the unsteady term with the reaction term, we obtain 


k (2 ) 


* * * * 

p C r dr 


HT 


h* 2 K* q / T (0.0) 

c c 


(2.46) 


Since both heat transfer with gas and catalytic chemical reaction 

can contribute to the heat-up process, one way of defining the total 

* 

heat-up time is 

1 = _J_ + _JL (2.47) 

* “ *( 1 ) * ( 2 ) 

T T T 

HT HT HT 

However, this formulation does complicate the nondimens ional 

equations. In the computations performed, we find that x^ is 

always smaller than x*' 2 ! Therefore, for simplicity, x ht is taken 

to be x* ll) as given by Eq. (2.45). With this choice. the 

HT 

dimensionless form of Eq. (2.44) becomes . 
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pC- 


a t 

d t 


rdr = r (T - T ) 

nk g c , s 


(2.48) 


+ 


* ' 

A 2 D 

c at 

* 

A (o,o) Nu 
g 03 



0.8 

y exp 

c 



( 1 - 



rdr 


where 


r 

nk 


Nu A 

x g 


Nu A (o,o) 

CO g 


(2.49) 


* * 

With r = t , we also get 

HT 


in gas phase 


6 = 
g 



« 1 


in catalytic layer 



1 , 



« 1 


*2 * 

h / a 

in substrate layer 5 = c s 

S * 

T 

HT 


Neglecting these terms with small parameter S’s, the dimensionless 
equations all become quasi-steady and are presented as follows : 
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Gas-phase 

continuity p u = p(o,o) uto,o) 


3 y 


-E 

*, >b , g 


species pu + J D ^ y g y c ,s^ ~ B i^ PY g ) exp( T * 


a t 


-E 

energy pu * — - + J (T - T ) = q B* (py ) b exp (-=2-) 

eneigy Q x H g c, s g 2 g 1 


equation of state 


P = p T 


(2.50) 


(2.51) 


(2.52) 


(2.53) 


Boundary conditions are given at the inlet of the reactor 


y (x=0) = 1 
g 

p (x=0) = 1 


T (x=0) = 1 (2.54) 

9 

u (x=0) = 1 


Solid-phase 

species in catalytic layer 


a 



D' 

am 


y° ' 8 exp ^ ( 1 - 



r 3 r 


0 


(2.55) 
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energy in catalytic layer 


3 

r 3 r 



) + D' 

at 


0 . 8 

y exp 

c 




= 0 


(2.56) 


The similarity between the energy Eq.(2.55) and the species 
Ecj . (2.56) suggests that a Damkdhler analogy (22] can be used to 
eliminate one of two equations in the computation. A detailed 
derivation of the Damkdhler analogy is carried out in Appendix B. 
The analogy is given by : 


T 

c 


T + 

c, s 


D 

at 


D 

am 


(y 


c, s 


y ) 


(2.57) 


Substituting Eq. (2.57) into Eq. (2.55) for T gives a modified 

C 

species equation as : 


a . a 
(r 3T 


r 3 r 


) = 


D y°‘ 8 exp|E 1- 

am c I c 


1 


D' 

r + — (y 

C » S r\/ C , S 


D' 

am 


y ) 

c 


(2.58: 


The energy equation for the substrate layer, Eq. (2.40), is reduced 
to 


3 


d_ 

a 


r 3 r 


(r 


) = 0 


(2.59) 
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Integrating Eq. (2.59) twice, the result is : 


T (r)= c Ln ( r ) + d 


where c and d are constants and it subjects the following boundary 
conditions : 


at r 


r ; 

c 


T = T 

s c 

a t a t 

^ s _ c 

X a r - a~r 


(2.61a) 


(2.61b) 


at r = r 

S 


a t 

s 

a r 


(2.61c) 


Applying Eqs . (2.61a) and (2.61c) gives 

T (r)= T (r ) (2.62) 

s c c 

that is, the substrate temperature is uniform in r-direction. Notice 
that it can be a function of time during transient heat-up . 

The energy-integral model represents an extension of th 
transient model in [5] where the catalytic layer is a surface (no 
structure). In the present transient model, finite catalytic layer 
thickness is included; temperature and species non-uniformity are 


allowed. 


CHAPTER III 


NUMERICAL SCHEME 

^ Full Transient Model Treatment 

For the gas-phase part which includes four Eqs. (2.34) to 

(2.37) and four unknowns — y^, T^, u, and p — the Euler explicit 
method with an upwind space difference is applied on the species 
equation and the energy equation for y and T ; the Euler backward 

<3 g 

implicit method is applied to solve the continuity equation for u. 
Finally, p is solved from the equation of state. Finite difference 
representations for the gas-phase equations are listed in the 
following: (superscript n+1 denotes quantity at new time level, and 

without n+1 denotes quantity at old time level n ) 


Gas phase 


species 


C = y-> + At Jdi 


gi-l J gl' S P (y c,si y q i 

<3 i 


(3.1) 


At * b- 1 h -F 

* — B , ?, y„ «xp ( — ) * y 

9 i 


gi 
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energy 


r 1 


t4i-) 

Ax 


A . J 

, At Hi f 'T T 1 ^ 

(T - T ) + -s (T - T ,) 

6 qi-l gi 5 p c,sl gi 

g <J 1 


u 


( 3 . 2 ) 


At * b - 1 b r E \ . f 

K p. y„, exp < — 1 ♦ T „ 


g S 2 r i ' g l 
g 


continuity 


n+1 


, n+1 x ( At x n + 1 n+1 l 

l«,< v p i 1 +< Tr )p . 


, _ n+1 n+1 x / At * 

(2p t - P 1 _ 1 H-r^) 


Ax 


( 3 . 3 ) 


equation of state 


n+1 


r,n+l 


( 3 . 4 ] 


Considering the Von Neuman stability analysis for the Euler explicit 
method, this method is stable when the following condition is 

satisfied : 


„ , At w m 

0 s (— — ) ( — - 
Ax * 


) s 1 


( 3 . 5 ) 
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Solid phase 

For the solid phase part, Eqs.(2.38) to (2.40) are solved by 
the implicit Crank-Nicolson method [23] which is an unconditional 
stable scheme with accuracy of 0[ (At ) 2 , (Ax) 2 ] A 
quasi-linearization technique (also called Newton linearization) 
[24] is implemented on the nonlinear heterogeneous chemical reaction 
rate to avoid numerical stability difficulties [25]. The resulting 
tridiagonal systems of linear algebraic equations which are 
diagonally dominant are usually solved by the LU-decomposit ions 
method [26]. The finite difference equations are derived in 
Appendix C and the following equations are obtained (where bracket 
denotes a matrix). 


species in catalytic layer 


M FI - M 

energy in catalytic layer 

MM • [",] 


(3.6) 


(3.7) 
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energy in substrate layer 

M [ C] - [ d 3 ] 

The order of computational execution for a monolithic reactor 
under cold start-up operation using the full transient model is 

(1) Specify the initial values of T g (x), y g (x=o), p(x), u(x) in the 
gas phase, and TU. r), y^x. r), T g (x, r) in the solid phase. 

(2) Solve the gas-phase equations (3.1) through (3.4) to obtain 

T (x) y (x), p(x), and u(x) at axial location Ax and new time 

g ’ g 

level n+1 . 

(3) Substitute the updated gas-phase properties from step (2) for 

A A a n D D Due to the 

matrix coefficients , A 2> A^, u 3 

application of quasi-linearization on the reaction rate, 
iterations are required in Eqs. (3.6) (3.7) (3.8) by using a 

simple iterative updating procedure. To do this, the profiles 
(r) T (r), T (r) used to start iterations are those of the 

previous time step, and the system solved for new y^r), T.(r), 

T (r) at the new time level At. This procedure can be repeated 

S 

iteratively until a convergence criterion is reached. 
(Agreement of successive iterations within a tolerance of 0.1 % 
of the dimensionless temperature is used in this study. ) 

(4) Through steps (1), (2). (3), the full transient model proceeds 

to the next axial location 2Ax at the new time level At. In 
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summary, computations are carried out by sweeping the x-axis 
from the inlet of reactor to the exit at every At, then 
advancing to next time step. 

Gas-phase equations (3.1) to (3.4) and solid-phase equations 
(3.6) to (3.8) are coupled, and they must be solved simultaneously. 
This complicates the problem considerably. However, the gas-phase 
properties are not sensitive to the variation in the solid phase, 
and a also small time step is used in computation. Therefore, the 
values of T" ^ and y"** in the gas-phase equations are approximated 
by T c,s and y c, s ’ respectively. The finite differences are 
evaluated on uniform grids in axial x direction and non-uniform grid 
in radial r direction. In x direction, it is divided by 31 points 

and gives Ax = 0.033. Because reaction occurs in the catalytic 

layer only, steep temperature and species gradient may exist in this 
layer but not in the substrate layer. Therefore, a smaller radial 
space step is used in the catalytic layer (Ar=0.05) and a larger 
grid step is used in the substrate layer. In the sample calculations 
to be presented, the catalytic layer is one third of the total solid 
thickness and the substrate is the other two third. Twenty one grid 
points are utilized in each of these layers. The time step At=0.1 is 
used in most cases. However, varied time steps are used for cases 
in which the gas ignition is strong. A smaller time step, At=0.01, 
is used during the gas ignition period. 



34 


3. 2 Energy Integral Model Treatment 

Rewrite the governing equations of the energy-integral model 

derived in chapter II. 

Gas phase 


continuity 


p u = p (0,0) u (0,0) 


(3.9) 


species 


a y 


p u 


— + J (y - y ) = B (py ) exp(-E /T ) 

x Dgc,s i g gg 


(3. 10) 


energy 


a T 


pu 


+ J (T - T ) 


6 X H g c,s 


= q B (py ) exp(-E /T ) 
n g 2 g g g 


(3.11) 


equation of state 


P = p T 


(3.12) 


B oundary conditions 

T (x=0) = 1 
g 


y (x=0 ) = 1 (3.13) 


p (x=0) = 1 


u (x=0 ) = 1 
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Solid phase 

The overall transient solid energy equation (2.48) is written in 
finite difference form : 






pC T ( t+At , r ) rdr = 


pC T(t, r) rdr + 


(3. 14) 


At 


{ 


r 

nk g 


T ) + 

c , s 


A 2 D 

c at 

* 

X (o, o)Nu 

g oo 


y c’ 8 e x p[ E c ( l- 4" ) ] rdr } 


The left hand side of the above equation represents the thermal 
energy In the solid (catalytic layer plus substrate layer) at time 
t+At. 


species in catalytic layer 


3 y 

3 -Cr c 


rdr Sr 


n' o. 8 

D y exp 

am c 


H 1 - 


T + 

c, s 


D' 

a t 

D' 

am 


(y - y T 

c , s c 


Boundary conditions 


(3. 15) 


at r = r ; 

d 


= r 


d y 

C 

d r 
a y 

c 

a r 


Sh (y - y ) 

c g 


(3. 16a) 


at r 


0 


(3. 16b) 
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Temperature profile in catalytic layer 


T = T + 

c c, s 


a t 


(y 


y ) 


c 


(3. 17) 


Temperature profile in substrate layer 

T (r)= T (r ) (3- 18 > 

s c c 

The overall transient solid energy equation (3.14) needs 
initial profiles for temperature T^r), T^fr), and species y c (r) at 
every Ax location along the channel. These profiles can be obtained 
from the full transient model results, assuming computation is 
switched from the full transient model to the energy-integral model 
at time t . 


Initial conditions 
catalytic layer 


substrate layer 


T (x, r, t ) 

c 0 

y (x, r, t ) 

c 0 


specified 

specified 


t 


) 

o 


T (x, r , t ) 

c c 0 


(3. 19a) 
(3. 19b) 


T (x, r, 


(3.20) 
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The computational procedures are listed as follows : 

(a) The overall transient solid energy equation (3.14) is 
integrated for one time step (At) at each x to find the new 


(b) 


(c) 


thermal energy integral 


pC T(t+At,r)rdr at that x location. 


Jr d 

The values at present time t are used in the right hand side 
of Eq. (3.14) in this calculation. 

To find the temperature profile which satisfies the energy 
integral at t+At, Eq. (3.15) and (3.17) are solved iteratively. 
The iterative procedure is as follows. First, estimating the 
boundary value T c g and Eq. (3.15) is solved by shooting 
technique [27] for y (r) subject to boundary conditions (3.16a, 
b). Once y c (r) is obtained, it is substituted into Eq. (3.17) 
for T c (r) ‘ The temperature profile is then used to evaluate 
the energy integral. If the value of the integral is not equal 
to that obtained in (a), a new T is estimated. This 

C, S 

procedure is repeated until convergence is obtained. The solid 
T s and y^ which satisfy the energy integral is the solid 
temperature and oxygen mass fraction at t+At. The process is 
repeated for all the grid points in the x-direction. 

With the surface temperature and oxygen mass fraction obtained, 
the gas-phase properties are to be updated to t+At. This 
process involved the solutions of quasi-steady equations (3.9) 
to (3.12). Since these equations are first order in x, the 
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Runge-Kutta method is used to integrate them from the channel 

entrance (x=0) to the exit (x=l). 

(d) The whole cycle is then repeated for the next time step. 


The shooting method reduces the solution of a boundary value 
problem to the iterative solution of an initial value problem. The 
usual approaches involve a trial-error procedure. In order to have a 
good initial guess, the first guessed y^ and is selected from 

the previous time step. The advantages of using this model are that 

(1) The Damkohler analogy eliminates the computation of the 
catalytic energy equation. 

(2) Assumption of " thermal ly-thin layer " eliminates the 

computation of the substrate energy equation. 

(3) Quasi-steady assumption also eliminates the computation of the 

continuity equation. 

(4) A large time step is allowed in computation. 


CHAPTER IV 


RESULTS AND DISCUSSIONS 


4.1 Introduction 

There are a large number of nondimensional parameters in the 
equations and in the boundary conditions derived in Chapter II. 
Associated with these nondimensional parameters are many time scales 
as shown in Table 2. Clearly the catalytic ignition problem can 
have a variety of complex behaviors depending on the choice of the 
values of these parameters. Instead of investigating the influence 
of all these parameters, in this study, we will limit our 
computational effort to fulfilling two basic objectives. 

The first objective is to understand the reactive structure 
variation with time in the catalytic layer during the ignition 
transient, since the new element of our model is the inclusion of 
catalytic layer in-depth reaction. The second objective is to test 
the energy-integral model as an efficient way to compute the heat-up 
transient of the igniter. 

To fulfill these two objectives, we need to be able to perform 
computation through the whole ignition transient using both the full 
transient model and the energy-integral model. Because of the 
smallness of the time step needed in the full transient method, the 
parameters of the test cases are deliberately chosen in such a way 
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that a reasonable computing resource is adequate. This means that 
the heat-up time of the solid is purposefully chosen to be small (as 
compared with values listed in Table 1). Nevertheless, the order 
between the various time scales is still preserved. 

4. 2 General Transient Behavior 

A representative case (referred to as Case 1) that illustrates 
the combustion processes from reactor start-up through transient to 
steady states are presented here. The values of properties for 
Case 1 are listed in Table 4 and the values of operating parameters 
are listed in Table 5. 

Figure 2 shows the axial profiles in the gas channel for gas 
temperature (T ), gas-phase oxygen fraction (y^), gas-phase reaction 
rate (W*), surface temperature (T ), surface oxygen fraction (y c g 

g c » s * 

) and maximum catalytic reaction rate (VM at times 10, 50, 100 x^ 

(where t* =500 x"). At t*= 0 the temperatures of the gas and the 
HT m 

solid are uniform and equal to the inlet temperature (400°K). The 

gas channel has been filled up with hydrogen and oxygen is 

* « * 

introduced at the inlet. At very early times (t =10x m =0.13x r ), the 
oxygen front has not been convected through the entire reactor and 
the oxygen concentration drops to zero abruptly at x=0.3 in the gas 
channel which can be seen in Fig. 2b. When the oxygen is flowing 
downstream, it is also absorbed by the catalytic layer. Oxygen is 
also consumed by the catalytic reaction which is represented by the 
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increased catalytic reaction rate in Fig. 2f. The gas temperature 
(Fig. 2a) is still at the level of the inlet temperature, but the 
solid temperature has a slight rise at the inlet (Fig. 2d) due to 
the catalytic reaction occurring there (Fig. 2f). Figure 3 shows 
the contour plot for T (x.r), T (x.r), y (x.r) and W (x.r) at time 

c s C C 

* * 

t =10 t^. Oxygen is absorbed uniformly in the radial direction and 
gradually propagates in the axial direction (Fig. 3b). Similar 
profiles in the solid reaction rate can be seen in Fig. 3c. 

Therefore, the non-uniformities such as species and reaction rate in 
the solid are primarily in the axial direction only. 

At slightly longer time (t =50 x ), the oxygen front just 

passes through the exit of reactor (Fig. 2b). The oxygen in the 
solid extends further in the axial direction, and non-uniform 
profiles appear near the inlet section (Fig. 4b). A concentrated 
hot spot can be seen in Fig. 4a in accordance with the local high 
reaction rate shown in Fig. 4c. As time increases (t =100 x 

m 

Fig. 2d), the solid temperature rises rapidly from a lower inlet 
temperature to a peak which is located slightly behind the entrance, 
then decreases monotonical ly until the exit of reactor. By 
contrast, the gas phase receives heat from the solid, and the gas 

temperature rises slowly and uniformly (t*=100 x*. Fig. 2a) because 

m 

of receiving heat from the solid phase and the convection effect. 

Because there is more oxygen in the gas phase near the tube 
entrance and also a larger heat and mass transfer coefficient, the 
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oxygen concentrations are higher near the entrance (Fig. 3b). It 
may be noticed that the solid temperature is non-uniform in the 
r-direction in the catalytic layer but is uniform in the substrate 
layer (e.g. Fig. 4b). The computed uniform substrate temperature 
from the full transient model is consistent with Eq. (2.59) because 
t *=100 x , it is much greater than the thermal diffusion time of 

m 

the substrate, and the substrate behaves as thermally thin. The 
reason that there is a temperature gradient in the catalytic layer 
is because of the existence of reaction. 

* 

As time increases, from one to two heat-up time , the 

* 

peaks in T (Fig. 6a) and T (Fig. 6d) grow larger. At time 2 
the gas-phase reaction becomes active near the exit of the reactor 
(Fig. 6c). As may be seen from Fig. 7 and 8, the catalytic reaction 
rate reaches the maximum value at time 2 x^ while the onset of 
gas-phase reaction occurs. It implies that the gas phase becomes 
active when the solid phase has reached maximum catalytic reaction 
rate. The oxygen has been consumed in most of the solid part except 
at the inlet region and near the surface in the downstream region 
(Fig. 7b and 8b). Likewise, the active solid region is shrinking to 
the inlet and the surface (Fig. 7c and 8c). The in-depth reaction 
zone in the inlet indicates that the reaction is kinetically 
controlled and the surface-confined reaction zone in the downstream 
indicates that reaction is diffusionally controlled. 

At steady state (t*= 3 t^), both the gas and the solid 
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temperatures reach the adiabatic flame temperature 1360°k (Fig. 6a 
and 6d) . The abrupt rise in T and W* (Fig. 7c) is due to the gas 
ignition. Once the gas is ignited, the gas temperature starts to 
exceed the solid temperature, and the heating of the solid is via 
energy transfer from the gas. This results in a slight temperature 
rise on the solid surface (Fig. 9a). Beyond the gas-phase ignition 
point (about two thirds of the bed length), the oxygen concentration 
drops to zero abruptly everywhere: the gas phase (Fig. 6b), the 
solid surface (Fig. 6e) and the in-depth layer (Fig. 9b). 
Therefore, the active catalytic reaction regions at steady state are 
further reduced to about two thirds of the total solid volume. The 
attainment of complete oxygen conversion in a short distance depends 
on the presence of homogeneous reaction in the gas channel; this 
special feature of catalytic combustor [28] also has been observed 
in experiments (in contrast to catalytic converter with only 
heterogeneous reaction). 

4.3 The Effect of Catalyst Loading 

Catalyst loading is defined as the weight of catalyst applied 
to the solid [29] . It is usually expressed either as a percent of 
total substrate weight or as weight per unit substrate volume. The 
influences of catalyst loading are twofold. First a larger catalyst 
loading leads to the faster heterogeneous reaction, and consequently 
the solid phase will generate more non-uniformities during the 
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transient period. Second, the catalytic ignition point moves toward 
the inlet as the catalyst loading increases. A reasonable catalyst 
loading of 1.8 kg/m 3 [30] has been selected in Case 1, and the 
computed results are shown in Fig. 2 to 9. In order to investigate 
the effect of catalyst loading, two catalyst loading values are 
assigned in Case 2 and Case 3 for 0.18 and 18 kg/m 3 , respectively, 
while the other parameters are the same as Case 1 (the catalyst 
loading 0* appears in Eq. (2.18)). Due to the lower catalyst loading 
in Case 2, the longer reactor length (20 cm) is needed to achieve 
both the catalytic ignition and the gas ignition in comparison with 

the shorter length used in Case 1 (5 cm). 

In Case 2, the transient responses in the reactor are very slow 

compared with Case 1, as shown in Fig. 10 at three time levels 5, 

7 15 x * The influences of reduced catalyst loading are obvious. 

’ ht' 

Since the initiation of catalytic reaction is delayed, the whole 

* 

combustion process in the reactor takes more time (15 t^) to reach 
steady state and a longer length (20 cm) to complete reactions. 

At t *_ 5 x * ( the solid and the gas temperatures are almost 
identical (Fig. 10a and lOd). Figure lOf and 11c show that the 
maximum catalytic reaction occurs at the exit but not at the inlet 
as in Case 1. The in-depth solid temperature (Fig. 11a) and species 
(Fig. lib) are uniform in the radial direction but not in the axial 

direction. 

As time increases to 7 t* the gas and the solid temperatures 

H 1 
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are rising at the same speed (Fig. 10a and lOd). Catalytic reaction 
rate increases and the peak region moves upstream as shown in 
Fig.lOf. The solid temperature is rising in the aft bed (Fig. 12a) 
and the solid species is decreasing similarly (Fig. 12b). 

At steady state (t =15 t ), the gas-phase reaction is 
triggered, which can be seen by a sudden rise in the gas temperature 
(Fig. 10a) and the gas reaction rate (Fig. 10c). Because of low 
reactivity, the oxygen penetrates all the way into the catalytic 
layer and there is little radial-direction gradient as shown in Fig. 
13b. The maximum catalytic reaction rate exists in the middle of 
reactor (Fig. 13c) and the gas-phase ignition point is located right 
after it. It implies that the solid phase and the gas phase 
response at the same speed. 

Case 3 has a higher catalyst loading which is ten times higher 
than Case 1 and one hundred times higher than Case 2 (see Table 3). 
The increase in the catalyst loading causes the catalytic ignition 
to occur at the inlet (Fig. 14f ) and also moves the gas— phase 
ignition point toward the inlet (Fig. 14c). Initially, there are 
strong radial gradients in the solid phase for T , y , and W* (Fig. 
15) due to a sudden catalytic ignition occurring at the front while 
the rest of the solid phase is not active yet. In addition, the 
catalytic reaction is so strong that oxygen is consumed near the 
solid surface and does not have a chance to diffuse into the solid 
(Fig. 15b). 
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As time increases (t =1 t ), the peak in T moves slightly 

HT c , s 

backward from the inlet (Fig. 14d), although the highest catalytic 
reaction rate is still at the inlet, and the active solid region is 
enlarged toward downstream (Fig. 14f). The solid temperature has 
less gradient in the radial direction (Fig. 16a) because most oxygen 
has been consumed by the surface layer (Fig. 16b). The maximum 
value of heterogeneous reaction rate is located at the inlet and is 
stabilized (Fig. 16c). 

At steady state (t*=3 x* ) , the gas ignition occurs in the 
middle of the reactor (Fig. 14c) , and the whole solid reaches the 
adiabatic flame temperature except in a very small inlet region 
where it has a slightly lower temperature (Fig. 14d). Oxygen in the 
gas phase is consumed quickly downstream of the gas-phase ignition 
point (Fig. 14b). Therefore, in the downstream region there is 
little reaction (Fig. 17c). The gas and the solid temperatures 
become equilibrized in this zone. 

Figure 18 gives the ignition boundary map as a function of 
reactor length and catalyst loading. Figure 18 shows that there are 
three different regions for ignition of monolithic catalytic 
combustor : a region with small reactor length and catalyst 
loading where no ignition occurs; a region with sufficiently large 
reactor length and loading factor where both catalytic and gas 
ignitions occur; and a region in-between where catalytic ignition 
occurs but not gas-phase ignition. Figure 18 shows that catalytic 
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ignition has to occur first before gas-phase ignition. It shows 
that at high loading, the catalytic ignition occurs at the tube 
inlet and the gas-phase ignition length becomes constant. But at 
low loading, the ignition length increases steeply with decreasing 
loading. For practical device (where there is a length limit), 
there will be a minimum catalyst loading below which the igniter 
does not function. Notice that Fig. 18 is plotted with a fixed 
incoming velocity and mixture ratio. In addition to its influence 
in the ignition boundary which is deduced from steady-state results, 
catalyst loading affects the transient response of the igniter. 

In summary, as the catalyst loading increases, the 
non-uniformities in the reactor, which include the gas phase and the 
solid phase, are diminishing at steady state. Concerning the time 
to reach steady state, if the catalyst loading is significantly 
reduced such as in Case 2, then it takes longer time. However, when 
the catalyst loading is increased, such as in Case 3 (0 =18 kg/m ), 

* * 3 

it still takes 3 t , as in case 1 (0 =1.8 kg/m ). This happens 
because the higher catalyst loading only strengthens a local 
heterogeneous reaction but the major parts of the solid are not 
affected. Therefore, the whole solid still takes the same time to 
reach steady state as a low catalyst loading case. 
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4 . 4 Th e Effect of Mixture Fuel/Oxidizer Equivalence Ratio 

The inlet equivalence ratio is a measure of the oxygen 
concentration in the mixture. In a fuel-rich mixture, when 0 is 
reduced toward unity (stoichiometric ratio), the heat of reaction 
per mole of mixture increases. This raises the adiabatic flame 
temperature which affects reactions in both the gas phase and the 
solid phase. 

Case 4 is given a higher oxygen concentration (0=4) than Case 1 

(0=8) and is used to investigate the effect of equivalence ratio. 

At time t*=0.2 x ( = 100 x ), the solid has a slight temperature 
HT m 

rise at the inlet (Fig. 19d), but the species on the solid surface 
are dropped suddenly (Fig. 19e) due to the strong heterogeneous 
reaction which can be seen in Fig. 19f. More detailed in-depth 
solid behavior can be obtained from the solid contour plot in Fig. 
20, The highest solid temperature is located at the inlet (Fig. 

20a) which is in accordance with the highest heterogeneous reaction 
rate shown in Fig. 20c. The in-depth species increases from the 
inlet to the exit (Fig. 20b). 

As time increases to t =1t » the solid ignites near the inlet 

HT 

r lx * Pie 19d), the gas phase is triggered after the solid 
HT ’ 

ignition position and a strong gas temperature rise can be seen in 
Fig. 19a. Two peaks appear in the solid surface temperature profile 
: the first peak is due to the catalytic ignition and the second one 

is due to the gas-phase ignition. After the gas phase is ignited, 
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the gas temperature exceeds the solid temperature after the gas 
ignition point. Therefore, a slight temperature jump appears in the 
solid surface temperature profile due to the heat transfer from the 
gas phase to the solid. However, before the gas ignition point, the 
heating process is reversed. Also these two hot spots on the solid 
surface can be seen clearly in Fig. 21a. As the gas ignition is 
achieved, the oxygen is consumed almost completely after the 
ignition point (Fig. 19b and 19e), and the heterogeneous reaction 
in the catalytic layer also ceases at the same location (Fig. 21c). 

By comparing the transient profiles in Case 1 and 3, we find 
that the effects of a more stoichiometric mixture are : (a) 
temperatures in the gas phase and the solid phase are higher; (b) 
less reactor length is needed to achieve the gas ignition; (c) 
there is less oxygen penetration into the catalytic layer; (d) less 
time to is needed to reach steady state (Case 4, 2 x* and Case 1 3 

HT * 

* 

T H1 .) • All those are due to the enhanced homogeneous and the 
heterogeneous reactions at higher flame temperature. 

4.5 The Cause of Non-uniformity in the Sol id 

In the model development, it has been assumed that the gas 
phase is one-dimensional and the solid is two-dimensional. In the 
above computations, we have discussed the cause of non-uniformities 
in the axial direction of the reactor. In this section, we discuss 
the non-uniformities in the radial direction of the solid. 
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Notice that associated with the solid energy equation and the 
species equation in the catalytic layer are the two Damkohler 
numbers Dat and Dam, and the other parameters such as q, Y^(o,o) , Sh 
and Be. First, concerning the species equation of the catalytic 
layer Eq.(2.22), the value of Dam is an important controlling factor 
which is the ratio between the mass diffusion time and the 
heterogeneous reaction time. Hence, when the mass diffusion time is 
slower than that of the mass consumption, a species gradient will 
exist in the radial direction. Similarly, in the energy equation of 
the catalytic layer, Eq.(2.25), the value of Dat, which is the ratio 
between the heat conduction time to the heterogeneous reaction time, 
controls the temperature gradients in the solid. When the heat 
conduction rate is slower than the heat release rate, temperature 
gradients will exist. In addition to conduction-reaction balance, 
the unsteady term is also a source for radial gradient. 

Case 3 is an example of increasing the catalytic reaction 
rate by increasing the catalyst loading ten times from Case 1. At 
early time (0.2 t* ) , a sequence of in-depth solid temperature 
profiles from the inlet (x=0) to the middle of bed (x—0.5) is shown 
in Fig. 23. The most non-uniform temperature profile exists at the 
inlet, and the peak is located in the middle of catalytic layer 
(r=0.3). Further downstream (x=0.005, 0.01, 0.5), the radial 
temperature profiles become more monotonic, decreasing from the 
surface toward the interior. Figure 24 shows the oxygen species 
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profiles, all monotonical ly decreasing from the surface as expected. 

The steeper gradient exists at the inlet. As time increases (1 x* 

ht’ 

Fig. 25), the hump-type temperature profile at the inlet disappears; 
a monotonical ly profile increasing from the surface toward the 
interior is evolved. These monotonic profiles indicate the approach 

toward a quasi-steady distribution. When the steady state is 

* 

reached at 3 r (Fig. 27), all temperature profiles are monotonic. 

Except at the inlet, the temperature distributions are almost 
uniform. The fact that the lowest temperature occurs on the surface 
of inlet is due to the higher heat transfer from the solid to the 
gas at that position. 

Another example is given by Case 4 in which the amount of 
oxygen concentration in the approaching mixture is double that in 


Case 1 . 

Consequently, the value 

of 

Da t 

increases 

due to 

the 

increase 

of heat of reaction q , 

and 

/ 

Dam 

decreases 

because 

of 


decreasing of Y^(o,o). This trend implies that the energy will be 
generated faster than the consumption of the species. In addition, 
the value of mass transfer coefficient Sh increases at the inlet 
region. Therefore, the species can diffuse into the solid more 

efficiently through the effect of increased Sherwood number. 

* * 

At small time (t =0.2 t ht ), a sequence of in-depth solid 
temperature profiles from the inlet (x=0) to the middle of bed 
(x=0 . 5 ) is shown in Fig. 29. At the inlet, the temperature profile 
is the opposite of the other locations (x=0.06, 0.12, 0.5). This 
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profile is due to the inlet has reached quasi-steady state. Oxygen 
is monotonical ly decreasing from the surface to the interior of the 
solid (Fig. 30). 

At larger time (t*=lx* ), the solid temperature is rising and 

HT 

the gradients in the radial direction are decreasing (Fig. 31). The 
inlet is the coolest part, and the rest of the solid has higher but 
uniform temperature. The temperature profiles reflect on the 
species profiles in Fig. 32, which shows a steep gradient at inlet 
and becomes small afterwards because of strong consumption. 

At steady state (t*=2x* ), the temperature profiles are similar 
as early time (t*=lx* ) but the temperature magnitudes are increased 
(Fig. 33). The species profiles are also similar as before and the 
magnitudes are decreased (Fig. 34). 


4.6 Computed Results Using Energy-Integral Model 

Although the full transient model is theoretically more 
accurate, it can be computationally costly. When the short 


transient period can be ignored, the energy integral 


method can 


offer an attractive alternative. There are two ways to apply the 


energy-integral model. One is to apply it right from the beginning 
of the transient, and the other is to first start with the full 


transient calculation and switch it to the energy-integral at a 


specified time. 

In this section, the computed results using the integral model 
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will be presented and compared with the results using the 
full-transient model. The purpose is to access the accuracy of the 
energy-integral model. In the first comparison, parameters of the 
Case 1 are computed using the energy-integral model initially. 

Figure 35 shows comparisons of the gas temperature at three time 

* * 

levels — 1, 2, and 3 t . At early time 1 t , the 

ht ' ht 

energy- integral model predicts lower gas temperature than the full 
transient model. However, at later times (e.g. 2, 3 x* ) , results 

HT 

from the two models compare favorably. Comparisons for the gas 

species are shown in Fig. 36. It indicates that the gas species is 
slightly over-predicted by the energy-integral model at 1 x* , but 

HT 

predictions from two models fall on one line afterwards. Hence, the 
predictions on the gas species profiles are reversed from the gas 
temperature. Figure 37 shows the comparisons of the gas reaction 

rates. Since the gas reaction becomes important after 1 x* , and at 

HT 

that time two models predict almost identical gas temperature and 
gas species, the comparisons of gas reaction rate have shown good 

agreement . 

To compare the detailed predictions in the solid phase, it is 
better to look simultaneously at the two dimensional contour 
plot in the solid and the in-depth plot for the radial profile. By 

comparing the solid contours in Fig. 7a with that in Fig. 38a at 1 

* 

x , We n °t* ce that although the oxygen and the reaction rate plots 
compare favorably, the temperature profiles predicted by two models 
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are somewhat different. Figure 39 for the in-depth solid 

temperature shows that the deviation is the smallest at the inlet 
(x=0) and increases afterwards. At the middle of bed (x=0.5), the 
slopes predicted by two models are different. This is because the 
aft bed reaches quasi-steady state later than the front part. (The 
energy-integral model implies a quasi-steady temperature 

distribution as shown by Eq. (2.57).) Thus, applying the 

0 Q 0 rgy— integra 1 model too early malces deviations at the aft bed 

section. Comparisons of the gas species profiles (Fig. 40) also 
show more deviation downstream, but the monotonic trend of y^ is the 
same , 

As time increases (t =2 x ), the catalytic reaction in the 

HT 

solid phase has reached maximum value. Therefore, the deviations 
from the two models’ predictions in T, y , W become smaller. The 

C C 

contour plots are compared by Fig. 8 and 41. For the solid 

temperature, both models predict a hot spot at inlet region; and the 
temperature distributions are very close. In addition, the in-depth 
temperature profiles compared in Fig. 42 show good agreement from 
the inlet to the middle of bed. Comparisons in the species are 

shown in Fig. 43; the two models predict identical profiles. These 

good agreements continue to exist up to the steady state as shown 
in Fig. 44 and 9. 

In summary, although the energy-integral model is applied at 

* 

the beginning, the comparisons such as T, y^, and W c have shown 
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great compatibility by using the simpler energy-integral model. The 
maximum deviation between two models predictions is the solid 
surface temperature at 1 x* — the value is 10 %. The maximum 

HT 

deviation in the gas temperature is 3 %. The saving in 

computational time has been achieved by using the energy-integral 

model in Case 1. Using Sun Spark Work station, Case 1 takes 52 

* 

minutes to reach steady state at 3 x in comparison to 30 minutes 

HT 

by the energy-integral model. If the heat-up time x* is increased 

HT 

100-fold, which is a more realistic value, the computational time 
using the energy-integral model should remain unchanged , but the 
computational time using the full transient model would be 5200 
minutes. 

In the next example, Case 3 is chosen to compare the results 
computed by the two models. The reason for choosing Case 3 is that 
it has a higher catalyst loading (/3*= 18 kg-Pt/m 3 ) which results in 
more non-uniform temperature profiles than Case 1 (/3*= 1.8 

3 

kg-Pt/m ). These constitute a more critical test for the 
energy-integral model which assumes quasi-steady state profiles for 
solid temperature and species; these profiles are monotonic. 
Moreover, due to the higher solid temperature gradient occurring at 
the inlet, more grid points are needed there. Hence, non-uniform 
and more grid points are applied in this case. 

Figure 47 shows the comparisons of the gas temperature. At 

« * 

early time (t =0.2 x ), the energy-integral model predicts lower 

HT 
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gas temperature in the middle of the bed, and both models predict 

consistently at the inlet and the exit region. As time increases 

(t*=lx* ), the energy-integral model predicts higher gas temperature 
HT 

* * 

than the full transient model. At steady state (t =31:^), 
comparison between two models is in good agreement. Figure 48 shows 
that the energy-integral model always predicts lower gas species, 
but the differences are small. Figure 49 shows the comparison of 
gas reaction rate. Both models predict the same maximum gas 
reaction rate and identical gas ignition position. 

* 

A comparison of the solid contour plot at 0.2 t ht in Fig. 50 
with those in Fig. 15 shows that the energy-integral model predicts 
higher solid temperature, lower species and narrower reactive 
region than the full transient model. These deviations can be seen 
more closely in Fig. 51. The relatively large discrepancy indicates 
that the profiles are still in the full transient region and that 
using the energy-integral model is not valid. 

As time increases to t =lx , deviations between two models are 
decreasing. From a comparison of Fig. 16 and 53, it seems that the 
energy-integral model predicts less temperature gradient in the 
radial direction. However, looking into the in-depth temperature 
profiles in Fig. 54, we can confirm that the deviations in the solid 
temperature are small. Near the inlet region (x— 0, 0.01), both 
models predictions fall on one line, but separate afterwards 
(x=0 . 5 ) . Comparisons of the solid species are shown in Fig. 55, and 
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are in good agreement. 

* * 

At steady state (t =3x ) , the comparisons of two models in the 

HT 

solid phase (Fig. 56 and 17) are similar as time t*=lx* . Because 

HT 

the solid phase has reached the maximum catalytic reaction rate 
after lx^. During 1 t* t to 3 t^., the aft part of the solid is 

receiving heat from the gas phase until the solid reaches steady 
state. These good agreements of the two models can be seen clearly 
in the solid in-depth profiles shown in Fig. 57 and 58. 

In summary, the comparisons of two models for Case 3 have 
shown the effect of catalyst loading on the deviations in the gas 
phase and the solid phase. The highest deviation in the gas 
temperature occurs at transient 1 x* , the difference is 4 %. The 

H T 

deviations in the solid phase are decreasing as time increases 

(maximum 15 % at 0.2 x ) and reach the minimum (less than 1 %) at 

steady state. Due to finer grid size (51 points) used in Case 3 
than Case 1 (31 points), computational time is longer in Case 3. 

The full transient models takes 74 minutes and the energy-integral 
model takes 46 minutes to reach steady state. 

In order to simulate a real monolithic reactor which has a 
longer solid heat-up time than those in the model’s cases, Case 5 is 
presented. In Case 5, the solid heat-up time is two seconds which 
is 100 times longer than that in Case 1. This is achieved by 
increasing the thermal inertia of the solid 100 times. The 

transient calculation is made by using the energy-integral model. 
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Figure 59 shows the computed transient results. Comparing the 
results in Fig. 59 with those in Fig. 6 (Case 1), we find that they 
are identical. These results are physically reasonable. They imply 
that the solid heat-up time is the correct time scale to 
nondimensionalize the transient equations. They also illustrate the 
advantage for using the energy-integral model. By using the 
energy-integral model, the computational time is identical in Case 1 
and 5. If the full transient model is used, the computational time 
needed in Case 5 will be 100 times longer than that used in Case 1. 

In order to have accurate transient results during whole 
combustion processes, the switching time from the full transient 
model to the energy-integral model should be greater than any time 
scales involved except the solid heat-up time. 


CHAPTER V 


CONCLUSION 

Two transient models for a monolithic catalytic igniter with 
in-depth heterogeneous reaction have been developed and solved 
numerically to study the ignition and transient phenomena. These 
models include a full transient model which accounts for the fastest 
unsteady event, such as mass diffusion across the catalytic layer, 
and an energy-integral model which accounts for only the slowest 
heat-up process in the transient. In this study, the reactants used 
are mixtures of rich-hydrogen/oxygen, and the catalyst is platinum. 
Two ignition phenomena are observed. First catalytic ignition 
occurs within the catalytic layer; later, gas ignition is triggered 
downstream in the gas channel. Once the onset of gas ignition 
begins, the reactor rapidly approaches steady state, and an abrupt 
gas temperature rise can be seen. Nevertheless, the solid 
temperature rises more uniformly due to the effect of low activation 
energy possessed by the catalyst and the large thermal inertia of 
the solid. During combustion processes, the front part of the 
reactor is the catalytic reaction controlled and the aft part is the 
gas reaction controlled. Moreover, initially the catalytic reaction 
is kinetically-controlled and becomes dif fusional ly-control led after 
the catalytic ignition is achieved. 
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Computation using selected model parameters shows that both the 
oxygen concentration and the temperature gradients can exist across 
the catalytic layer. The oxygen concentration always decreases from 
the surface toward the interior. But temperature profile can be 
nonmonotonic in the reactive layer during the initial fast 
transient. At a later time, the temperature distribution behaves in 
a quasi-steady manner and becomes monotonic. 

The effect of catalyst loading has been investigated, the study 
indicates a large effect on the ignition length and ignition delay 
time. An ignition map using reactor length and catalyst loading has 
been constructed; the map shows the existence of three regions : no 
ignition, catalytic ignition only, and catalytic and gas-phase 
ignitions. 

Comparisons between the calculated results using the full 
transient model and the energy-integral model show that after the 
initial short transient period (greater than all the time scales 
except the heat-up time), the two models produce similar results. 
Therefore, unless short transient information is needed, the use of 
the energy— integral model will yield correct results and save a lot 
of computation time. If initial transient is important, the most 
logical approach is first to use the full transient model and then 
switch to the energy-integral model. 
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APPENDIX A 


DETERMINATION OF CHEMICAL REACTION RATES 


A. 1 Introduction 

The combustion processes in a monolithic reactor include two 
parts: one is a homogeneous reaction taking place in the bulk gas 

phase; another is a heterogeneous reaction occurring on the catalyst 
which is deposited on the solid wall. In this study, 

r i ch-hydr ogen/oxygen mixtures are used as reactants; the catalyst is 
platinum. In order to simplify the model formulation, one-step 
global reaction rates in the gas and the solid phases are desirable. 

Concerning the homogeneous chemical reaction occurring in the 

gas phase, the reaction rate is deduced from the concept of ignition 

delay. The parameters of the reaction rate are calibrated by 

experimental data. Experimental data also have shown thatthe 

reaction rate depends on deficient oxygen only. On the other 

hand, for the heterogeneous chemical reaction occurring in the solid 

phase, several experimental studies have been reviewed. These 

studies give a consistent activation energy 5500 cal/mole, and the 

. , . 0.8 . 

reaction rate depends on the oxygen concentration P q only. 
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A. 2 Gas Phase Reaction Rate 

We propose a one-step global model for rich-hydrogen/oxygen 
reaction. The proposed model is based on the spontaneous ignition 
theory [10], and the parameters of the global rate are calibrated by 
the experimental data. The overall energy equation for an adiabatic 
reaction within the gas channel is written as follows. 


* 


P c 


♦ d T 
P d x 




* * 

= q W 


(A. 1) 


A one-step global reaction rate is assumed as. 


* » * a * b » m 

W = K P P q exp (-E /R T ) [kg-0 2 /m 3 /sec ] (A. 2) 


* 

Let T be the initial temperature of the mixture and assume that the 

O 

preignition reactant consumption is negligible and T*- T*« T*. 

o o 

Nondimensional parameters are : 



* * 


R T 



(A. 3) 


O 
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One approximation is made for exponential term in reaction rate as 
follows : 


e xp ( — — ) = exp 
R T 


[ 1 

L R (T -T + T ) J 


(A. 41 


exp 


# * 

* T - T 

-L. (1 + _°) 


R T 


1 j * exp ^ 0(1-7)) j 


Substitute Eq.(A.4) into Eq. (A.2), and integrate Eq.(A.l) 
after arrangement : 


1 



e q W (T ) 

( - ° - ) t 


* * 

P C 

p 


(A. 5) 


0-p 

When the mixture is ignited, i.e. e => °°, the critical time is 
called the ignition delay time t . Therefore ^ is defined as: 


* * * * 

p C R T 

p 2_ 

• • • * 

E q W (T ) 

O 


(A. 6) 


The one-step global reaction rate Eq.(A.2) has four parameters: 

„ * 

rate constant K, exponents a and b, and activation energy E . They 

will be determined by fitting the computed ignition delay time into 
the experimental data! 11, 12] which supply the ignition delay time 
for mixtures with various equivalence ratio, pressure and 
temperature. Therefore, the reaction rates from the fuel-rich to the 
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stoichiometric case can be calibrated. Procedures for calibration 
are discussed in the following section. 


* 

Activation energy E 

Differentiation of Eq.(A.6) with the neglect of T* relative to 

O 

* * 

E /R provides a way of obtaining the activation energy. 


d (In t ) 

i 

d (1/T* ) 



(A. 7) 


where E is the slope of plot t , vs. 
directly from the experimental curve. 


1/T* 

O 


and value 


* 

E 


can be read 


Exponent of oxygen b 

Choose two mixtures (denoted by 1 and 2) with different 
equivalence ratios, and their ignition delay times are measured 
under the same temperature and pressure. Then the exponent of 
oxygen b can be deduced from the following equation : 


t 

ii 


* 

t 

12 


Cp p 

1 ^ ol 


Cp. 


o2 


(A. 8) 
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Exponent of pressure a 

Choosing one mixture, two ignition delay times are measured at 
the same temperature but different pressure level. Then the 
exponent of pressure a can be defined. 


* 

t 

11 


* 

t 

12 



(A. 9) 


Rate constant K 

* * * 

Substitute these three parameters E , a , b into Eq.(A.6) for 

a specified mixture for which the ignition delay time is known, the 

* 

value for rate constant K can be determined. 

* 

T 

Criti cal temperature c r 

Unfortunately, one global reaction rate can not fit the 
experimental data over a wide temperature range; the reason is due 
to the chain reactions involving H0 2 [13] having opposite effects on 
the ignition only at low temperature and high pressure conditions. 
Therefore, two global reaction rates are required to fit the whole 
temperature range, and the switching point is determined by the 
critical temperature T* . For specified pressure and oxygen 
concentration conditions, if we equalize these two global reaction 
rates, a critical temperature T can be obtained. 
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Table A. 1 lists the corresponding parameters for two sets of 
gas-phase reaction rates. Fig. A1 shows that the computed ignition 
delay time based on two derived reaction rates is in good agreement 
with the measured ones in both fuel-rich and diluted stoichiometric 
cases . 

A. 3 Catalytic Reaction Rate 

A number of experimental studies [14, 15, 16] to determine the 
overall rate of rich-hydrogen/oxygen reaction occurring in a porous 
catalytic layer have been reviewed. Even though they used a 
different amount of catalyst loading and size of catalyst, and 
operated at a different temperature range, the conclusions are the 
same :(1) the catalyst reaction rate depends on the concentration of 
oxygen only; (2) the exponent of oxygen concentration is 0.8 and the 
activation energy is 5500 cal/mole. The reason to choose the 

selected reaction rate reduced from Maymo, et. al [14] is because it 

* 

is based upon the catalyst loading 0 which is defined as the total 
amount of catalyst deposited on the solid. Therefore, it will give 
more flexibility to study the influence of thecatalytic reaction 
rate on ignition phenomena. 
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A. 4 Summary 

A set of chemical reaction rates, the unit given as 
[kg-0 2 /m 3 /sec] , has been derived, and it will be used in this study. 

Two reaction rates are required in the gas phase, and the switch 

* 

point is determined by a critical temperature T^, which is the 
interception point between these two rates. In this study we have 
T* =1200 °K. 

c r 


gas-phase reaction rates 


when T ^ T 

cr 


1 8 - 1 . 5 # 0. 1 m m 

V/* = 4 x 10 P p exp (-79480/R T ) 


when T ^ T 

cr 


* * * * * 

W = 9 x 10 P p exp (-19780/R T ) 

H o 


(A. 10) 


(A. 11) 


catalytic reaction rate 


W* = (3*1.772 (p* T*) exp (-5230/R T ) 

c o 


where (3 is the catalyst loading [kg-pt/m 3 ] . 


(A. 12) 


APPENDIX B 


DAMKOHLER ANALOGY 


Rewrite the steady state equation (2.55) and (2.56) in the 
catalytic layer which have been derived in chapter II. 


species in catalytic layer 


o 


a f g y c 

r d r r dr 


D y°* 8 exp(-E /T ) 

am c c c 


(B. 1) 


energy in catalytic layer 


0 = 


r d r 


(r 


d T 

c 

d r 


) + D y°' 8 exp(-E /T ) 

at c c c 


[B. 2 ) 


The similarity between Eq. (B.l) and (B.2) results in the following 
equation. 


-j-j 

r d r| 



am 


T + D* 

c at 


y 

c 



(B.3) 
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Integrate Eq. (B.3) once : 

r\ * / 

r (D T + D y ) = constant (B.4) 

d r am c at c 

Integrate Eq.(B.4) from the catalytic layer surface r d to the 
in-depth distance r. 


d (D T + D y ) = constant 

am c at c 


(B. 5 ] 


D T 

am c 


+ D y 

at c 

r 

d d 


= 0 


(B.6] 


Equation (B.6) is subject to the following boundary conditions : 


at r = r ; T = T 

d c c, s 


(B. 7a] 


y = y 
c c, s 


(B. 7b] 


Solving Eq . (B.6) results in the Damkohler Analogy 


T = T + 

c c» s 


a t 


( y - y ) 

c , s c 


(B . 8 ] 


am 


APPENDIX C 


DERIVATION OF FINITE DIFFERENCE EQUATIONS 


In the full transient model, the solid-phase equations which 
include species, energy in the catalytic layer and energy in the 
substrate layer are solved by the implicit Crank-Nicolson scheme. 
The finite difference equations are derived in the following steps. 


species in catalytic layer 
transient species equation (2.38) 


T 


g 


p y 

O C 

~d~r 


a 

r d r 


(r 



r / °- 8 t ' 

- D y exp E 

am c c 



(C. 1) 


Each term of Eq.(C.l) is discretized by using the central 
difference. The subscripts i and j denote the axial and the radial 
grids , respectively. 


p y 

o c 


n+ 1 
'cj • 


cj 


a t 


At 


(c. 2; 


a y c 
a r 2 


24r 2 t 


n+1 ^ , n+1 

y + y - 2( y + 

cj+l cj+l cj 


x n + 1 

y )+ y + y 
c] cj-1 7 cj 




(C . 3) 
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a y c 
a r 


n+1 n+1 

/ - y + y + y 

c j + 1 cj-1 cj + 1 cj-1 

4 r Ar 
j 


(C.4) 


A quasi-linearization technique is implemented on the nonlinear part 
of catalytic reaction rate, which is obtained from the Taylor series 
with second and higher orders neglected. By using the form 


f (T"*\ y"*‘)= f(T, y) * [14 • ~H ] 


T n+1 - T 


n+1 

y - y 


(C. 5 ) 


The nonlinear part of reaction rate is assigned as 


-E 

_ . 0.8 , c . 

f (T , y ) = y exp( ) 

c c c 1 


(c.6: 


The linearized reacion rate can be obtained 


f (T n + 1 , y" + 1 )= exp( A >[ y' 8 (-2 - ) + (— 

c c lie i 

c L c i 


y‘ 8 E 

c c j ^n+1 
2 c 


, 0 .2 v n+1 

+ ( . 8 y ) y 

C C 


(C. 7 ) 
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Boundary condition at r = ( where denotes j=l 


9 y 


cl 


a r = Sh ly c r V 


(C. 8 ) 


The second order derivative term at r = r is obtained from the 

d 

Taylor series 


a y 3 y .2 

y = y + ( — 3 -^)Ar + ( ^1) + 0 (Ar 2 ) 


c 2 cl a r 


a r 


2 2 


(C. 9 ] 


a 2 y 


cl , 


Solving Eq.(C.9) for ( ), we obtain 

3 r 2 


a 2 y 


cl 2 


8 r 2 Ar 2 


y - y (1+Ar Sh) + Ar Sh y 

c2 c 1 g 1 


(C. 10) 


Similarly, at r = r , the boundary condition is (where denotes j=N ) 


a y 


cN 


a r 


= 0 


(C. 11 


a 2 y 


cN 2 


2 .2 

a r Ar 


[ y c«-l - y <» ] 


(C. 12) 


Substituting Eq.(C.2) through Eq.(C.12) into Eq. (C.l) for j = 1, 2, 
• • • -N, the following N simultaneous algebraic equations can be 


obtained . 
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where 


where 




with 




(C. 13) 


[A^] represents the following tridiagonal matrix: 



(C. 14) 


y n+1 and D represent the following vectors 
C 1 J 


r] - 


i — n+l 
y cl 
n+l 
y c2 


n+l 

L- y cN 


and 


M- 


u 


12 


IN 


I c. is; 


O' = 



(C. 16) 


D 

am J 


At 


/ 

D exp E fl- 
am c 



(C. 17) 
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-1 A r * A c* \ AtSh 

i = -=r- “ 4<r(l+ArSh) + 

ill r 

g i d 


.8 D y°’ 2 

ami cl 


(C. 18) 


b = 4cr 
N 


(C. 19) 


“1 " 0 2 

a = - 4cr - . 8 D y 

IN T amN cN 

gi 


(C. 20 ) 


nH + 1 


D = -=- -(4crAr- )Shy + 

n T r gi 

gi d 


D «iy;![- 2 --r (1 --f } 

L cl cl J 


(C.21) 


D = -^ N + D y’ 8 

IN T amN cN 

gi 


[■*- 


E T n+1 

~Y (1- ~Y±-) 

cN cN 


] 


(C . 22 ) 


when j - 2, 3, . . . N-l 


a 


ij 



gi 


.8 D 

am J 


0. 2 
'cj 


(C. 23) 


r = r + ( j-1 ) Ar 
J 


b 

j 


= cr 



At 


4Ar 


c 

j 


1 

= cr + 

r 

j 


At 


4Ar 


r ^ r s r (C.24) 

s j c 

(C . 25 ) 
(C. 26 ) 


D = -b y + 2cr 
ij J J-i I 


T 

gi 


+D 

a m j 


(. 


2 - 


- )y * 2 ly - c s 

cj J^cj j J 


c J + 1 


(C. 27 ) 


+ D 

a m j 
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energy in catalytic layer 

Applying the same technique on the energy equation in the 
catalytic layer, the following N simultaneous algebraic equations 
can be obtained. 




A 

2 




(C . 28 ) 






a 

c 

0 



21 

l 




b 

a 

c 


-1 


2 

22 

2 


A 

= 

0 

b 

a 


2 



3 

23 







c 





N-l 




. b 

a 


— 


N 

2N — 1 


r T n+1 n 

cl 

^n+1 

A c2 

and 

D 

2 

— 

D 

21 

D 

22 

T n+1 

— cN — 




D 

2N 


( C . 29 ) 


(C . 30 ) 


with 
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at j 


At D exp E (1" 


a t 


L 



(C.31) 
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a = 
21 


-6 - 


c 


4cr (l+ArBc) + 


AtBc 


r 

d 


+ 


D 

a tl 


Y '! E 

cl c 


T 


2 

cl 


a = 
2N 



4 <x 


A (4cr 



At 

Ar r 

i 


-) 

c 


// 


+ D 

atN 



E 

c 


T 


2 

cN 


£ 

D = -5 T - ( 4<rAr + — )BcT - d” y' 8 [ .2- — - ■ 

21 c cl r gi atl cl T 

d L c 1 


n+ 1 

.8 y cl " 


cl J 


At 


D = -5 T -(4crAr + 

2N c cN r 


) A 


r n+1 

s 2 

"a F 


- D y 

atN 


■•[. 2 - 

cN[ 


cN 


n+ 1 

. 8 y cN I 


cN 


when j = 2, 3, . . . N-l 


2 j 


-5 - 2cr + 

c 


at j 


y' 8 e 

c J c 


C J 


D = 
2 J 


b T 

C j- 



T 

J 


•C T 

C j + 


D 

atj 





n+1 
. 8 y c j 

y cj 


where Ar^ and Ar^ are the radial space steps used in the 
and the substrate layers, respectively. 


( C . 32 ) 


(C. 33) 


(C. 34 ) 


(C. 35) 


(C. 36) 


(C. 37) 


catalytic 
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energy in substrate layer 


[ *3 ] [ C] ’ [ “3 ] 


(C . 38 ) 
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APPENDIX D 


COMPUTER PROGRAMS 


D. 1 Flow Chart 
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D.2 PROGRAMS 


C 

c 

c 


c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROCRAM TRANS 


CATALYTIC COMBUSTOR START-UP TRANSIENT 


REAL Lf.MMK.MOLO 

C0MM0N/THER1/ YHO(31).YO(3l).YH(3l).YOSai),YHOV(31), YH0SV(31). 

♦ YlC(3l).YLl(3l).XJH(31),RNK(3l).WG(3l).RO(3l).T5(3l).T(3l). 

♦ XX(3l).DXt3l).DY.DT.R0I.YHI.Y0t.TI.CP,DS.XLS,XLDS.RT(3l) 
C0MM0N/THER2/B(21).G2I).Bl(21),G(21).B3(21),C3(21).RR(21).RS(21) 
C0MM0N/THER3/SH(3I).BG31),YRR(3L21).TSI(3l.21),DAM.DAT 

♦ , Y0SU3l.21).TS2(3l.21),WS(31.21).DEU3l) 

COMMON/THER*/ TgNEW(3l ) . Y0NEW<3l).R0NEW(31). UNEW( 31) . Uf 31 ) 
common/th erS/wsmix(3U.tsrtn( 31) 


OPEN! UNIT-9 .fTLE-'OUr.STATUS-'UNKNOWN*) 
OPEN! UN IT-10. FILE »'OU2’.STATUS-‘UNKNOWN’) 
0PEN(UNrT-ll.FlL£»0U3‘. STATUS- ‘UNKNOWN’) 
OPEN! UN rr- 12. TILE- ’ 0U4‘ .STATUS-'UNKNOWN’ ) 
OPEN! UNIT-13. nL£-'0U5‘ .STATUS-'UNKNOWN* ) 
OPEN! UNIT- 1*. FILE- ‘ 0U6‘ .STATUS-'UNKNOWN’ ) 
OPENtUNIT-lS.nL£-‘OU7’ .STATUS- ‘UNKNOWN*) 
OPEN! UNIT- 16.FILE* ‘ 0U8’ .STATUS-‘UNKNOWN* ) 
0PEN(UNTT-17.nLE-‘0U9‘. STATUS- ’UNKNOWN*) 
OPEN! UNIT-18. FILE- 010’ .STATUS-‘UNKNOWN* ) 
OPEN! UNIT- 19. FILE- ' 011‘ . STATUS- 'UNKNOWN* ) 
OPEN! UNIT-30. FILE- • datl‘ ) 


DATA PF.LF.3S.N .NY / 209.6 . 1100 . , LE-a. 31 . 21 / 

DATA XLEW1.W0.WH.CC / L . 32. . 2. . .08205 / 

DATA XPR.AOAS.xki / .719 .2. . .0* / 

DATA CL.AL.BL.EL.DKRO / 4.E18 . -L5 , 0.1 . 79*80.. S7800. / 

DATA CH.AH.BH.EH.DELT / 9.E5 , L . L . 19870., .99/ 

DATA R00T1.R00T2.X0UNT/ 0. . L , L / 

DATA ALPH1.CC2.ES /l.E-5 . 2.5 . S320. / 

DATA VOID.FRA / .5 , .67/ 


INPUT PARAMETERS 

REACTION 1. A H2 * 02 « 2 H20 * { A-2 ) H2 (CAS PHASE) 

REACTION 2. 2H2 * 02 - 2H20 (SOLID PHASE) 

EL.EH.ES - ACTTV. ENERCY FOR CAS (LOW. HICH T). SOLID (CAlVMOL) 
XLS . DS - BED LENGTH (Ml. CHANNEL HYDRA UL. OLA. (M) 

DEF - EFFECTIVE DIFFUSIYTTY OF 02 IN CATALYST (M2/S) 

SCC . SC - HEAT CAPACITY OF CATALYST. SUBSTRATE (CAL/KG/K) 

SRC . SR - DENSITY OF CATALYST. SUBSTRATE (KG/M3) 

BVt - LOADING FACTOR 

DY. FRA - RADIAL CRID LENGTH. RATIO OF CATALYTIC LAYER 
THICKNESS / HALF WALL THICXNESS 

DT. DTS - TIME STEP IN FULL TRANSIENT. ENERCY INTECRAL MODEL 
KMAX. KST - MAX COMPUTATION TIME. TIME TO SWITCH MODEL 


READ 130, 91 SRC. SRC. SCC. SX. SR. SC DEF 

READ (30. 9) DS, XLS. RM. XPI. TI. UIUS. BV1 
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READ (30.10) DT. DTS. FRA. DY. KST. KMAX 
9 FORMAT(7E10.3) 

10 DORMAT(4E10.3.2U0) 

C 

c OPERATING PARAMETERS AT INITIAL UPSTREAM (I.U.) 

C TI . UI » I.U. TEMP. IK) . VEL. (M/S) 

C XPI . ROI » I.U. PRESSURE (ATM.) . DENSITY (KG/M3) 

C CP GC » I.U. SPECI. HEAT (CAL/KG/K) . GAS CONST. (M3- ATM/KG/S) 

C YHI . YOI - I.U. MASS FRACTION OF H2 .02 

C XKI . XPHI » I.U. CONDUCT. (CAL/M/S/K) . HZ/02 EQUIV. RATIO 

C ALPH1 * 02 IN CATALYST 

C XLEWl. - LEW. NO OF 02 



UI » UIUS*(1.0*A0AS)/A0AS 

c 

C •••• THERMO. DATA 
c 

CALL THERMO (RM.DHRO.MOLO.MMIXTAB.DHR.AA.XPHI.CC.XPI) 



C ”*• BREAK TEMP FOR GAS PHASE TWO-EQ MODEL 

CALL BREAK (CL.AL.BL.EL.CH.AH.BH.EH.TCR.MOLO.DHR.XPI) 

C ***• GAS PHASE REACTION CONSTANTS 

CL » CL/XPI**(-AL)*(ROI*YOI)”BL 

EL » EI^TI/l. 987 

EH - EH/TI/1.987 

QV - DHR/CP/TI 

GR1- XLS/ROI/UI 

XLE123- XLEWf.667 

XMU - 1.05E-5*(TI-373.)*l.SE-8 

REY « UIUS*DS*ROI/XMU 

ALPHA * XKI/ROI/CP 

AD2UL * ( ALPHA/D S/DS)/(UI/XLS) 

C **** SOLID PHASE CONSTANTS 

E2 - ES/TI/1.987 
QV2 » DHR/SCC/TI 
THIC*.433*DS/A0AS/2. 

RL1 ■THIC'FRA 
RL2 -THIC’(l.-FRA) 

RI1 »0S/2./RLl 

RI2 -RI1*L 

RIO *RI2*RL2/RL1 

DAT -BV*QV2*RL1**2. ‘SCC/SKC 

DAM »BV*RL1**2./DEF/Y01/R0I 

C 

C •*** MATRIX COEFF TO SOLVE TS L YOS EQ 



CALL COEFF (N.NY.NR.NT.RI1.RI2.RI3.DY2.DTY.X2.X4.D2.DT2.22.Z4.SKC. 

♦SK.RK) 


ENERGY INTEGRAL EQ. CONSTANTS 


AP —DAT/DAM 

H 2 *0TS*DAT*SKC*2-/XXI/3. 66 

RC *<SR*SC*<RI3**2. -R12**2. ) ♦ SRC*SCC*( RI2**2. -Rn**2. ))/2. 
R0C1«SRC*SCC/RC 

R0C2*SR*SC*(RI3 # *2. -RI2**2. J/2./RC 

C 

C *••* TIME SCALES 

C 

tres*xis/ui 

TMC »V0ID*RL1**2./DEF 
TTC *RL1 **2. *SRC*SCC/SKC 
TTS *RL2**2- *SR # SC/SK 
tts2*rll # *2. # sr*sc/sk 
XLE2-TTS2/TMC 

DIFF*2. •V0ID*RLl**2./XLE123/3. 66/ ALPHA 
T0TAL*2. *RC*RLi**2./3.66/XKI 

C •••• INITAL CONDITIONS (CAS £ SOLID PHASE) 

CALL IN IT (N.NY.NR, GC.MMIX, UI,TMC, GCC.DTX, DLT, CR3, DMASS, TIME) 

C 

C — INPUT PRINTOUT 

C 

CALL WINP (RM.XPLUIUS,KMAX,SKC,SRC,SCC,DEF.BVCC,KP,SK.SR,SC,RLL 

+RL2.TCR,XPHI.MOLO f TAB.R£Y,TMC.TTC.TTS,XL£,XL£2,TOTAL.N,NY,NR, 

+0 IFF, NT) 

C 

C TIME INTEGRATION 

C 

XNUO » 885*SQRTtA0AS*XPR*REY) 

S55 CONTINUE 

ROOT * 5MR00T1 ♦ R00T2) 
lF(R00T.L£.(U)WU-.0444/R00TO.46-L34*R00T 
IF(ROOT. GT.O.i)YNU»O.Oli/ROOTO. 66 
YROOT»YNU-XNUO 

IF (ABSIYROOTLLT.O.Ol) CO TO 666 
IF (YROOT.GT.O.O) ROOTi-ROOT 
IF (YROOT.LE.O.O) ROOT2*ROOT 
CO TO S55 
666 CONTINUE 

KOUNT*KOUNT*l 

KOU-KOUNT-l 

if (kou .1L fast) then 

TIME-TIME+OT’TMC 

else 

1 1 me»ume-^lts *total 
endif 

C 

C •— CAS PHASE EQ. (FULL TRANSIENT EQ. AND EXPLICIT METHOD ) 

C 

if (kou .It. kst) then 

CALL CAS (N.QV,CRI.CL,BUEL,CH,BH.EH.TCR,XLEI23.CCC.XPLDTX.DLT. 

* ROOT, AD2ULJQCI ) 
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tts2*rtl ## 2. *zr* sc/sk 
XLE2-TTS2/TMC 

DIFF-2. *V0ID*RLl , *2-/XL£123/3. 66/ ALPHA 
TOTAL-2. *RC*Rn~2./a 66/XXI 


** INITAL CONDITIONS (CAS l SOLID PHASE) 


CALL INrr (N,NY,NR.GCJi4MIX t UI.TMC,GCC,DTX,DLT t GR3,DMASS,TIME) 


— INPUT PRINTOUT 


CALL WINP (RM,XP!,UIUS,KMAX.SKC,SRC,SCC,DEF,3VCC,KP,SK,SR,SC.RL1, 
♦RL2,TCR,XPHLMOLO,TAB.REY,TMC.TTC,TTS,XLE,XLE2,TOTAL.N,NY,NR. 
♦DIFT.NT) 


•••■ TIME INTEGRATION 


XNUO » 385 B SQRT(AOAS*XPR*REY) 

555 CONTINUE 

ROOT * 5VROOT1 ♦ ROOT2) 

IF( ROOT. LE.0.1) YNU- . 0444/R0OT+3. *6- L 34 • ROOT 
IF( ROOT. CT. 0. 1 ) YNU*0. Oll/ROOT +3. 66 
YROOT* YNU-XNUO 

IF (ABS(YROOT).LT.O.Ol) GO TO 666 
IF ( YROOT.GT.O.O) ROOT! -ROOT 
IF (YROOT.LE.O.Q) ROOT2-ROOT 
GO TO S55 
666 CONTINUE 

KOUNT-KOUNT+i 
KQU-KOUNT-l 
if (kou .It, kst) then 
TIME -TIME HTTHMC 
else 

time*tinie+dts "total 
endif 

—* GAS PHASE EQ. (FULL TRANSIENT EQ. AND EXPLICIT METHOD ) 
if (kou ,1 l kstJ then 

CALL GAS (N,0V,GR1,CL,BL,EL,CH*BH.EH,TCR,XLEI23,GCC,XPLDTX,DLT. 
* ROOT, AD2UL, XXI ) 


GAS PHASE EO. (QUASI-STEADY EQ. AND RUNCE -KUTTA METHOD) 


else 

CALL QUASI {N t QV.DMASS,GRLCL,BL,EL,CH,BH,EH,TCR,XPHLXLEi23. 
♦ KOU.XPLROOT,AD2UL,XKI) 
endif 

GAS PHASE REACTION RAte 


DO 700 I*t,N 
ELL-EL/TU) 

EHH*EH/TII) 

IF {(Til) .LE. TCR ). AND. (ELL.LT. 100. )) THEN 
WC(I)-CL*(RO(I)*YO(I))* # BL/EXPfEL/Tn)) 

ELSE IF {(TCI) .LE. TCR LAND. (ELL. GE. 100. )) THEN 
WG(I)*0. 

ENDIF 

IF ((TCI) -CT. TCR LAND. (EHH.LT. 100.)) THEN 
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SUBROUTINE THERMO (RM.DHRO.MOLO.MMIX.TAB.DHR.AA.XPHI.GC.XPI) 

C 

C~ THERMO DATA 

C 

REAL MMIX.MOLQ 

COMMON/THERI/ YH0{31).Y0(31),YH(3UY0S(3U, YHOV(3l).YHOSV(3L). 

♦ YIC(31),YLI(3l),XJH(3U,RNK{3l),WC(3lLRO(31LTS(3U,T(31) t 

♦ XX(3i).DX(3l),DY,DT.ROI,YHI,YOLTI.CP,DS,XLS,XLDS.RT(3l) 

C •** DHR * HEAT OF REACTION (CAL/KC-02J 

C CP * SPECIFIC HEAT OF MIXTURE (CAL/KG-K) 

C TAB - ADABATIC TEMP (K) 

C AA - NO. OF MOLES FOR H2 

C XPHI * EQUIVALENCE RATIO 

C Al.BL.Cl - COEFF. OF CP FOR H2 
C A2.32.C2 - COEFF. OF CP FOR H20 

C A3.B3.C3 - COEFF. OF CP FOR 02 

C 

DATA A1,81.C1,A2,B2.C2,T0/ 6.947, LE-4, 1.6E-7, 7.256, 1.149E-3. 

♦ 9.43E-3, 298./ 

DATA A3.B3.C3/ 6.148. 1.S51E-3. 3.07E-7/ 

C 

F(X)*(AA-2. )*{A1*(X-TO)-B1*(X*X-TO*TO)+CI*(X* # 3.-TO # *3.)H’2.*(A2* 
+<X-T0)+B2*(X*X-T0*T0)+C2*(X**3. -T0**3. ))-CC 
C 

AA*16. *RM 
MOLO*l./(AA*U 
MMIX*MQLQ # 32. >( I. -MOLO ) *2. 

YHI*2./MMIX*(1. -MOLO) 

YOl*YHl/RM 

XPHI*RM*S. 

ROI*XPI*MMIX/CC/TI 

XLDS*XLS/DS 

C CP AND DHR 

CPOl * 6. 148-0. lE-3*TQ-9. 23E-7*T0*T0 
CPH1 * 6.947-2.0E-4TO^i.3lE-7TOTO 
CPHOl * 7.256*2.3E-3*TO+2.S3E-7TO*TO 
CP02 * 6.148>3.1E-3TI-9.23E-7TITI 
CPH2 - 6.947-2.0E-4Tl + 4.8iE-7TI , TI 
CPH02 - 7.256>2.3E-3*TI*2.33E-7*Tl*TI 

dhh - Ai*nr-T0)-Bi*rnTi-T0*T0)-Hr4*nTT3.-T0 # ru 
DHHO * A2Tn-T0)^B2*nrn-T0T0I^C2Tn # r3. -T0 # T3. ) 

DHO - A3*m-T0)+B3*(TI*n-T0*T0)-C3*(n**3. -T0**3. ) 

DHR * (DHRO+OHH+.5*DHO-DHHO)*LE3*2./32. 

CP - (M0L0*CPO2+(i. -MOLD )*CPH2) # 1000. /MMIX 

C 

C ADIABATIC FLAME TEMP CALCULATION 

CC-2. •DHRO^AA'DHH^DHO 
XI* 100. 

XF *10000. 

A-FCXI) 

9 B-F(XF) 

IF (A*B) 99.99.40 
99 WW*.5*(XF+XI) 

IF (ABSIWW-XIM.E-3) 40.40,41 
4L YI-FIWW) 

IF (YI-A) 30.30,31 
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30 XF-WW 
CO TO 99 

31 XI-WW 
CO TO 99 

40 TAB-WW/TI 
RETURN 
END 


SUBROUTINE BREAK (CL.AL.BL.EL.CH.AH.3H.EH.TCR.M0L0.DHR.XPl) 

C 

C THE BREAK POINT FOR THE CAS PHASE REACION RATE 



REAL MOLO 

COMMON/THER1/ YHO(31).YO(3l).YH(3l).YOS(3l).YHOV(3l).YHOSV(31), 

* YICI3U. YLI(3l).XJH(3l),RNK[31),WG(31),RO(31).TS(31).T(31). 

* XX(3l),DX(31),DY,DT,R0I,YHI,Y0I.TI,C7,DS,XLS.XLDS,RT{3l) 

C 

F(X)*C1+C2*AL0C(C3*X)+C4/X 

C 

CI-ALOG(CL*ELyCH/EH*XPI**(AL-AH)*(MOLO*32.)**(BL-BH)) 

C2-8H-BL 
C3-< . 0820S/XP!) 

C4-(EH-EL 1/1.987 
XI-lOO. 

XF-l.ES 

A-F(XI) 

9 B-F(XF) 

IF (A*8) 99.99.40 
99 WW-.S^XF-XI) 

IF (A8S(WW-XI)-l.E-3 ) 40 . 40.41 
41 YI-F(WW) 

IF (Y1*A) 30.30.31 

30 XF-WW 
CO TO 99 

31 XI-WW 
CO TO 99 

40 TCR-WW/TI 
RETURN 
END 
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SUBROUTINE COEFF (N.NY,NR,NT.R!1.RI2.RD.DYZ.DTY.X2.X4.DZ.DTZ.22, 
♦ Z4.SKC.SK, RK) 



C •••• CONSTANTS FOR THE SOLID PHASE 


COMMON/THERi/ YHO(31).YO(31).YH(31).YOS(3l).YHOV<31).YHOSV(3l). 

* YIC( 31 ), YLI{3l),XJH(31).RNK(31).WC(3l).RO(3l),TS(31),T(3I), 

♦ XX(3l).DX(31),DY,DT,ROLYHI,YOLTl.CP.DS.XLS.XLDS.RT(3l) 
COMMON/THER2/B(21),C(21).Bl(21),CM21).B3(21),C3(21).RR(21).RS(21) 
XLLN*DT/2./DY/DY 

BETT-DT/4./DY 
DTY*OT/DY 
X2 *2. *XLLN 
X4 »4. *XLLN 
DO 100 J»l.NY 
RT(J)*n.OAT(J-l) , DY 
100 RRU!*RT(J)*DY*RI1 
DO 200 J-2.NY-1 
B(J) »XLLN-BETT/RR(J) 

C(J) »XLLN*BETT/RR(J) 

BUJ) »8(J) 

200 CKJ) -C(J) 

Ctl) «-X4 
CU1) — X4 
3(NY) »-X4 
31(NY)— X4 

C 

C SUBSTRATE LAYER CONSTANTS 


DZ*2. *DY 

NR»i RI3-RI2)/DZ*l 

IF (NR .EQ. I) THEN 

02*R13-RI2 

NR-2 

ENDIF 

XLLN2=0T/2./D2/D2 
BETT2*0T/4./DZ 
DTZ-OT/DZ 
DYZ-OY/DZ 
22 -2. ‘XLLN2 
24 «4. *XLLN2 
NT-NY-NR-l 
DO 300 J-l.NR 

300 RSU)»FL0AT(J-l)*DZ*RI2 
DO 400 K»2.NR-1 
B3(K) -XLLN2-BETT2/RS(K) 

400 O(K) *XLLN2<‘BETT2/RS(K) 
DO SOO J-l.NR-l 

SOO RT(NY*J)«FL0ATU)*DZ*1. 
C3(l) — Z4 
B3(NR)— Z4 
RK-SKC/SK 
RETURN 
END 
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K&SSST* 

irsssss 

C0*MM0N/THER3/SH(3l).BC(311.YRR(31.21),TSl(31.21).DAM,DAT 

*Co2N/ : S^ ( Tg^(3!).YONEW(31).ROI^W(31).UNEW(3l).U(31) 

real mmix 


initial 


CONDITIONS FOR QUASI-STEADY CAS PHASE EQUATION 


DO 100 1*1, N 
T (I)*l. 
u (I)*l. 

100 R0(I)*1. 

YOU) *1. 

YH(1) -l. 

YH0(1)*0. 

R0(1) «l. 

DMASS-R0(1)*U(1) _ 

C^Nm^C^mON FOR FULL TRANSIENT CAS PHASE EQU ATION 


TIME * 0. 

Y0NEWUJ-1. 

TgNEW(I) *1. 

RONEWUl-t. 

UNEW (l)*l. 

gcc*cc/mmix 

dlt*xls/ui/tmc 

0TX*0T/D0X 

CR3*0TX/DLT 

C^NmAL CONOmONS OF THE CATALYTIC ^NO S UBSTRATE LAYER 


DO 200 1-l.N 
Tsm *i. 
Yosm-o. 
DEL(I)-.9S 
DO 300 J-I.NY 
TSl II.JM. 

300 Y0S1(I.J)*0. 

DO *00 K*l.NR 
400 TS2 (I.K)-l. 
200 CONTINUE 
RETURN 
END 


o n n 
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SUBROUTINE WINP (RM.XPI.UIUS.KMAX.SKC.SRC.SCC.DEF.BVCC.KP.SK.SR.SC 

* .RL1.RU.TCR,TICN,XPHI,M0L0,TAB,R£Y.TMC.TTC.TTS,XLE.XLE2 TOTAL N 

♦ NY.NR.DIFF.NT) 


" INPUT PRINTOUT 


COMMON/THER1/ YHO(3l).YO(31). YH(3l), Y0S(3U.YH0V(3l).YH0SV(3l). 

* YIC(31J,YLI(3U.XJH(3l),RNK(3U. WC(3I),RO(3I).TS(3l).T(3l). 

* XX(31).DX(3U.DY.DT.ROI,YHI,YOI.TI,CP,DS,XLS.XLDS.RT(3I) 
COMMON/THER3/SH(31).BC(31),YRR(3l,21).TSl(3I.21),DAM.DAT 

* , Y0Sl(3l.21),TS2(3l,21). WS(3l,21).dei(3l) 

WRITE(I9.12) 

WRITE(19.22)RM.XPHI.XPI.UIUS.ROI.KP,KMAX 

WRITE(19.13) 

13 F0RMAT(/3X. TT , 3X. 'TAB' .7X. 'TICN' . 6X. ’TCR' ) 

WRITE! 19. 33 m.TAB.TICN.TCR 
WRITE119.23) 

23 FORMAT(/3X,'TMC',7X. 'TTC',7X. "TTS’,7X.’DIFF’,6X, ’TOTAL’) 

WRITE! 19.33 ITMC. TTC, TTS, DIFF, TOTAL 
WRITE! 19. 15) 

IS FORMAT(/3X.'K-S*,7X,’RO-S’.6X,'CP-S\6X.’DAM\7X.’DAT\7X.'XL£' 
*,7X.’XL£2’) 

WR ITE! 19. 33 1SK. SR. SC. DAM, DAT. XLE.XLE2 
WRITE! 19. 17) 

17 FORMAT(/3X. ’REY’ . 7X. ’MOLO’ . 6X. ’ CP’ . 3X. ’DT* ,3X, ’ DS’ . 3X. ’XLS’ ) 

WRITE! 19. 33 1REY, MOLO , CP.DT. DS, XLS 
WRITEU9.14) 

14 FORMAT!/ 3X.’K-C’.7X.’RO-C,6X.’CP-C.6X.’DEF’.7X.’RU’.7X.’RL2’) 
WRITE!19.33)SKC, SRC.SCC.DEF.RL1.RU 

WRnE!19.18) 

18 F0RMAT(/3X. ’N* , 9X, 'NY' . SX. ‘NR’ . 8x. ’NT ) 

WRITE! 19.19 IN.NY.NR.NT 

19 F0RMAT14110/) 

22 FORMAT (5E10.3.4X.I6.4X.I6) 

33 FORMAT (7E10.3) 

RETURN 
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SUBROUTINE CAS (N,QV.GR1,CL,BL,EL.CH,BH.EH.TCR.XLE123,GCC.XPL 

* DTX.DLT.ROOT.AD2UL.XXI) 

- FULL TRANSIENT GAS-PHASE EQUATION SOLVED BY EULER UPWIND EXPLICIT 
METHOD 

COMMON/THERI/ YHO(31).YO(3n,YH<3l).YOS(3l).YHOV(3l).YHOSV<31). 

* YICI3U. YLI(31).XJH(31),RNK(3n, WC(31),RO(3I).TS(31).T(3l). 

* XX( 3 I).DX( 3 U.DY,DT,ROI.YHI.YOI,n i CP.DS,XLS.XLDS,RT( 31 ) 

COMMON/THER4/ TgNEWDl). YONEW(31),RONEW(31).UNEW(3l).U(3I) 


— CORRELATION OF HEAT i MASS TRANSFER 


DO 100 t»l.N 

XK=«( . 04* (T( I ) TI/273. ) ‘ . 94 l/XKI 

XST»XX(I)*2.*AD2UL*XK-ROOT 

IF (XST .LE. .1) XNU».0444/XST*3.46-1.34*XST 

IF (XST .CT. .1) XNU-.011/XST *3.66 

RNK(I)-XNU*XK/3.66 

XJH(I)*4. *XNU*.AD2UL*XK 

IF (I .EQ. 1) CO TO 999 


•* CAS PHASE EQ -- T (X). YO (X). RO (X). U (X) 


c 


CR2 a DLT*R0(I) 

ELL*EL/T(I) 

EHH*EH/T(I) 

IF ((TCI) -LX- TGI LAND.lELL.LT.lOOH THEN 
Wl*CU*CRl # (RO(I)*YO(I)) # *9L/EXP(El^T(I)) 

ELSEIF (CTtl) -LX. TCR LAND. ( ELL. GE. 100)) THEN 
Wl-O. 

ENDIF 

IF t(T(I) .GT. TGI LAND.IEHH. LT. 100)) THEN 
Wl»CH*ORl - (RO(I} - YQ(I))**BH/EXP(EH/T(I)) 

ELScIF ((T(I) -GT. TCR LAND.(EHH.GE.100)) THEN 
Wl*0. 

ENDIF 

TGNEW(I)*T(I)HJ(I)*DTX/DLTTni-l)-T(I)WXJH(I)*(TS(I)-T(I))> 

>QV*Wl)*DT/GR2 

XJD-XJH(I)*XL£123 

Y0NEW(I)«Y0(I) + U(I)*DTX/DLT*(Y0(I-1)-Y0{I))^(XJD # (Y0S(I)-Y0(I))- 

>Wl/Y0I)*DT/GR2 

IF (YONEW(I) -LT. 0.) Y0NEW(I)»O. 

RONEW(I)»XPl/GCC/ROl/TI/TGNEWfn 

UNEW( I )*(DLT*(RO( I )-RONEW( I ) }+OTX*UNEW( I- 1 ) *RONEW( I) )/DTX/( 2. * 
^RONEW(I)-RONEW{I-i)) 

999 CONTINUE 
100 CONTINUE 
do 200 i*t,n 


tU)»tGnew(i) 

yo(i)*yonew(i) 

ra(i)*ronewti) 


200 u(i) «unew(i) 
RETURN 
END 
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SUBROUTINE FULLT (NX,NY,NR.NT.E.RI1.R12.K0U.DTY,X2.X4,KMAX.Z2.Z4, 
* DTZ.DYZ.DZ.XLE.XLE2.KP, RK.SKC.DEF.SK. TIME.TOTAL,8V. tab) 


full transient equation of species and energy eq. 

UNIFORM GRID SYSTEM ( CATALYTIC l SUBSTRATE LAYER ) 


DIMENSION A(21).A1(21).YSS<21).YS2(21),DE121).EM<21).EM1(21>.TNEW 

* (2I).YNEW(21),YX(21).YXl(2l).A3(21).EM3<21),TNEWZ(21).YX3(21). 

* XL(21).D(21).U(21).TSMAX(3l) 

COMMON/THER1/ YHO(3l), YO(3l).YH(3l).YOS<31).YHOV(31).YHOSV(31). 

* YIC(31).YLI(31).XJH(31),RNK(3l).WC(3l).RO(3l).TS(3l).Tl3l). 

* XX131) DX(3U.DY.DT,ROI.YHI. YOI.TI.CP.DS.XLS. XLDS.RTI31) 
COMMON/THER2/B(21),C12l).Bl(21),Cl(2l).B3(2l).C3(21).RR(2l).RS(21) 
COMMON/THER3/SH(31),8C(3l).YRR(3L2I).TSU31. 21 J.DAM.DAT 

* , YOSU3l.21).TS2l31,21), WS(3l,21),DEL(31) 
common/ther5/wsmax( 31 ). lsrm(31 ) 

DO 1000 M.NX 
DO 100 J-l.NY 
TNEWU1-TS1 (I.J) 

YNEW(J)-YOSlU.J) 

IF (YOSUI.J) -LE. 0.) THEN 
YS2(J)-0. 

YS8UI-0. 


ELSE 

YS2(J)-l./Y0Sl(I.J)**-2 

YS8(J)”Y0SUI.J)"*-3 

ENDIF 

100 OE (J)»OT/EXP( E/TSUI. J)) 
DO 200 J - LNR 
200 TNEW2 (J)» TS2 (I.J) 


C CONSTANTS FOR CATALYTIC LAYER 


DO 300 J-2.NY-1 

A(J) — l./T(t)-X2-.3*DAM"YS2(J)"DEUI 
300 All J)— ■■XL£*X2»DAT*YS3( J)*DE( J)*E/TS1(LJ)**2- 

A ll)-l./T(I)*X4*U.*OY*SH(I))-OT*SHU)/RIl*OAM*DEm*.8*YS2U) 
A(NY)»X4*0AM*DE(NY)* 3*YS2(NYM./TII) 

A1(1) M XL£*X4"(1.»BCII)"DY)-DAT"DEU) , YS8(1)*E/TS1(I.11**2.-BCU)/ 


-RI1*DT 

Al(NY)— 0AT*DEtNY)*YS8(NY)*E/TSl(I.NY)**2.*XLE-DTZ/RK/RI2+X4 , ( 

* l.*OYZ/RK> 




c CONSTANTS FOR SUBSTRATE LAYER 


DO 400 J-2.NR-1 
400 A3U) — XLE2-Z2 

A3(l) ■XL£2*Z4"(l.*RK/DY2)-RK"Cmf/RI2 

A3(NR)-Z4*XL£2 

c .... FINITE DIFF. METHOD FOR CATALYTIC LAYER (SPECIES EQ.) 

ITER-0 
ITER1-0 
666 CONTINUE 

IF ((ITER .EQ. 20) -OR. UTERI -EQ. 20)1 THEN 
WRITE!*.*)’ KOU'.KOU.’ I (CAT) ** M 


ENDIF 

DO SOO J-2.NY-1 

EM(J)-“3(J)*Y0S1(I. J-1)*Y0S1(I.J)*(X2-1./Tn)+0AM*YS2IJ)*0E(J) 


*(.2 
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♦-£/TSl(I,J)))-Y0St(I.J*l) , CU)+DE(J)*DAM*E*YS8(J)/TSl(I.J)**2.» 

500 ♦TNEW(J) 

EM(l)»Y0Sl(I,l)/T(I)+( X4*DY-DT/RIl)*SH(I)*Y0(I)-DAM*DE(l)*YS8(I)* 

♦ (.2-E/TSl(I.l) # (l.-TNEW(l)/TSl(I.I))) 

E.M(NY)»-DAM*DE(NY)*YS8(NY)*(.2-E/TS1(I,NY)*(1. -TNEW(NY)/TS1(I NY) 
♦)HYOSl(I.NYVTtI) 

CALL MATRLU ( A . B . C , XL , 0 . U . NY) 

CALL SOLVER (XL , D . U . EM . YX .NY) 



C •••• FINITE DIFF. METHOD FOR CATALYTIC LAYER (ENERGY EQ.) 

DO 600 J-2.NY-1 

EM1( J)a-Bl(J)*TSllI.J-I)+TSl(I. J)*(X2-XLE)-TSI(I.J*1)*C1(J)-DE! J)* 

600 •►DAT*(YS8(J)*!.2-E/TS1(I.J))*.8*YX( J)*YS2(J1) 

EMl(I)»TSl(I,l) , XLE*(X4*DY-DT/RIl) # 8C(I)*T(I)+DAT*DE(l)*(YS8(l)* 

♦( . 2- E/TSH 1. 1 ) )♦. 3* YX( 1 ) • YS2( I ) ) 

EMl(NY)*DAT"DE(NY)*(YS3(NY)*(.2-E/TSl(I,NY))*.8*YX(NY)*YS2(NY)) 
XLETSU I , NY )*TNE W2( 2 ) *( X4 *0 Y+0T/RI2 1/RIODZ 
CALL MATRLU (AI . Bl. Cl, XL , D . U . NY) 

CALL SOLVER (XL , D . U . EMI . YX1 . NY) 

C 

ERR1-0. 

ERR2*0. 

DO 700 K-I.NY 

ER1 -ABS( YXOC)-YNEW(K)) 

ER2 -ABS(YXHK)-TNEWOO) 

ERR1-AMAXKER1.ERRI) 

ERR2-AMAX1(ER2.ERR2) 

YNEW(K)«YX (K) 

700 TNEW(K)-YXHK) 

IF ((ERR1.GT. l.E-3).0R.(ERR2.GT. LE-3)) THEN 
(TER»fTER*l 
CO TO 666 
END IF 



C •*** FINITE DIFF. METHOD FOR SUBSTRATE LAYER (ENERCY EO.) 

DO 800 J-2.NR-I 

800 EM3(J) »-33(J)*TS2(U-I)+TS2(I.J)*(22-XLE2)-TS2(I.>l)*C3(J) 

EM3( 1 ) »TS2( 1.1 ) *XL£2+< Z4*D2-DT/RI2 ) *RK/DYTNEW(NY-1) 
EM3(NR)»XL£2*TS2(I.NR) 

CALL MATRLU (A3 . B3. C3. XL . D . U . NR) 

CALL SOLVER (XL . 0 . U , EM3 . YX3 . NR) 

DO 900 K-l.NR 
900 TNEW2(K)«YX3(K) 

C 

C •••• CHECX TEMP. GRADIENTS AT THE INTERFACE 

C 

ER1*ABS(TNEW(NY)-TNEW2( 1 ) ) 

IF (ER1.GT. S.E-2) THEN 

ITERI-ITER1+1 

ITER «0 

GO TO 666 

ENDIF 

999 CONTINUE 

C 

C MAX SOLID REACTION RATE 

C 

wmax*0. 

DO 910 J*l t NY 
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WS(U)« BV*YNEW( J)**.8/EXP(E/TNEW( J)) 

REF2- WS (I.J) 

WMAX- AMAX1 (REF2.WMAX) 

TS1 (I.J)- TNEW(J) 

910 YOSlU.J)- YNEW(J) 
wsmax(i)*wmax 

C 

C — MAX SOLID TEMP IN RADIAL POSITION 

C 

call maxx (ny,tnew,ian) 

TSRM(I)»RT{KM) 

tsmax(i}-tnew(lun) 

C 

DO 920 J-l.NR 
920 TS2 (I.J)- TNEW2(J) 

TS (I)-TSl (1.1) 

YOSd)-YOSl(U) 

1000 CONTINUE 

C 

C— END OF FINITE DIFFERENCE 

C 

C — - MAX GAS and solid TEMP IN AXIAL POSITION 

C 

call maxx (nx,t,L) 
call maxx (nx,tsmax,Li) 

C 

C — MAX SOLID and gas reaction rate IN AXIAL POSITION 

C 

call maxx (nx,wsmax,m) 
call maxx (nx,wg,mi) 

RETURN 

END 


subroutine maxx (n,yy,LL) 

c*“** to find the max. value 

dimension yy(31) 
gre-yy(l) 

LL-l 

do 100 i-l t n-l 

if (jry(i+D .gt. gre) then 

gre-yyU+l) 

U«(+l 
end if 

100 continue 
return 
end 



non 
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SUBROUTINE MATRLU (A.B.C.XL.D.U.N) 


— MATRIX SOLVER FOR LU-METHOD 


DIMENSION A(21),B<21).C(21),XL(21).D(21).U(21) 
D(1)»A(I) 

uuj-cu) 

DO 20 I-2.N-1 
XLII)*B(I)/D(I-I) 

D(I)*A(I)-XL(I)*U(I-l) 

UU)-C(I) 

20 CONTINUE 

XL(N)»B(N)/D(N-l) 

D(N)»A(N)-XL(N)*U(N-1) 

RETURN 

END 


SUBROUTINE SOLVER (XL.D.U.Q.X.N) 
DIMENSION XL(21).Q(21).D(21).U(21).X(21).Y(21) 
Y(1)*QI1) 

DO 20 I-2.N 

20 Y(I)«Q(I)-XL(I)*Y(I-1) 

X(N)-Y(N)/D(N) 

IF (X(N) -LT. l.E-5) X(N)«0. 

DO 30 J«N-1.1.-1 
X(JMYU)-UU)*X(J*I))/D(J) 

IF (X(J) LT. I.E-S) X(J)-0. 

30 CONTINUE 
RETURN 
END 
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SUBROUTINE QUASI (N.Q,DMASS,GR1.CL.BL.EL.CH.BH.EH,TCR,XPHI,XLE123, 
+ KOU.XPI.ROOT.AD2UL.XXI) 

C 

C— QUASI-STEADY CAS-PHASE EQUATION 

DIMENSION XXY(4).XKZ(4),XXT(4) 

COMMON/TH ERL/ YHO(31).YO(3l).YH(3l).YOS<3aYHOV(31).YHOSV(3I). 

* YIC(3l), YLI(3l),XJH(31),RNK(31),WG(3l).RO(3l),TSl3n.T(3I). 

* XX(31),DX(31),DY,DT,ROI, YHI,YOl.TI,CF.DS.XLS,XLDS,RT(3t) 

C — CORRELATION OF HEAT l MASS TRANSFER 



DO 100 I-l.N 

XK»( .04*(T(I)*TI/273. !**.94)/XKI 
XST-XXm*2.*AD2UL*XK*ROOT 
IF (XST ,L£. .1) XNU*. 0444 /XSTO. 46-1. 34 *XST 
IF (XST .GT. .1) XNU=.01l/XST 0.66 
RNK(I)«XNU*XK/3.66 
XJH(I)»4. *XNU*AD2UL*XX 
C 

IF (I .EQ. N) GO TO 999 
M«l 

DO 120 L»t.4 

ELL-EL/TT 

EHH*EH/TT 

IF (TT .LE. TCR ) THEN 
W1«CL*GR1*(R0(I)*Y02)**BL/EXP<EL'TT) 

ENDIF 

IF (TT .GT. TCR ) THEN 

W1»CH*CR1*(R0(I)*Y02)**9H/EXP(EH/TTI 

ENDIF 

XKY( L )-( -XJH( I ) •XLE123*( Y02- YOSI I))- Wl/YOD/DMASS 
XKZ(L)-( -XJH( I)‘XLE223*(YH02-YH0SV(I))*Wl/YHI*36./32. J/DMASS 
XKTt LM -XJH( I) *(TT-TS( I ) )*Q*Wl)/DMASS 
120 CONTINUE 

YO(I-I) -YO(I) *DX(I)"(XKY(l)«-2. , XKY(2!*2.*XKY(3)*XJCY(4))/6. 

IF (YO(I*l) .LT. 0.) YOII+D-O. 

YH0(I*l)-YH0(I)*0X(I) , (XXZ(l)*2.KK2(2)*2.*XKZ(3kXlCZ(4»)/6. 

IF (YHO(I*l) .LT. 0.) YH0IMM3. 

T (Ifl) *T(I) ♦OX(I) , (XKTU)*2.*XrT(2)*2.*XKT(3)*X]CTI4))/6. 

IF (T(I*l) .LT. L) T(I-l)-l. 

YH( MR -( L Y0( I+U l/XPHI 
100 RO(I*l)-l./T(I+l) 

999 CONTINUE 
RETURN 
END 
C 

SUBROUTINE ROUTIN !I.M. YH02. Y02.TT.XKY.XKZ.XKT) 

DIMENSION XKY(4).XXZ(4).XKTI4) 

C0MM0N/THER1/ YHO(3U.YO(3I).YH(3l).YOS(31).YHOV(3U,YHOSV(3U. 

* YIC(3l).YLI(3i).XJH(3l).RNK(31).WG(3l).RO<31>.TSI31).Tl3n. 

* XX(3I).DX(31).DY.DT.ROI.YHI,YOI.TI.CP.DS.XLS.XLDS,RT(31) 

IF (M .NE. I) GO TO 102 

YH02«YH0(I) 

Y02 -YO(l) 

TT -Ttl) 

M »M*l 
GO TO 119 
102 CONTINUE 

IF (M .NE. 2) GO TO 103 
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yhoz-yho(i)*dx(i)»xkz<i)/2. 

Y02 =YO(I) ♦ DX(I)*XKY(l)/2. 
TT »T(I) +DX(I)*XKT(l)/2. 

M =M+1 
GO TO 119 

103 CONTINUE 

IF (M .NE. 3) GO TO 104 
YH02 s *YH0(I)+DX(I) , XKZ(2)/2. 
Y02 »YO(I) +DX(I)*XKY(2)/2. 
TT -T(I) ♦DX(I)*XKT(2)/2. 

M -M+l 
GO TO 119 

104 CONTINUE 
YH02=YH0(I)*DX(n , XKZ(3) 
Y02 -YO(I) ♦DX(I)*XKY(3) 

TT =T(I) +DX(I)*XKT(3) 

119 CONTINUE 

IF (YH02 .LE. 0.) YH02-0. 

IF (Y02 .LE. 0.) Y02 =0. 

IF (TT .LE. 0.) TT -0. 

RETURN 

END 
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SUBROUTINE SHOOT! NX.NY,NR,NT,E.H2.RI1.RI2,K0U,R0CI.R0C2. UMAX, KP, 

* BV.AP.KST.DTS.total.tiroe) 

C TO SOLVE ENEXCY INTECRAL EQ. BY SHOOTING METHOD 
C NONUNIFORM YOS ,TS , NON-UNIFORM GRID 

C 

DIMENSION TTSUM<3l).TSR(21).YSR(21).YT(21).CJ(21).csmax(3l) 
COMMON/THER1/ YHO(31),YO<3l),YH(3I).YOS(3l).YHOV(3l). YHOSV(3l). 

♦ YIC(31),YLI(3l).XJH(31),RNKl31),WG(3U.RO(31),TS(31),T(3U. 

* XX(3l),DX(31),DY,DT,ROI,YHI,YOI,TI,CP.DS.XLS,XLDS,RT!3l) 
COMMON/THER3/SH(3l),BG3l).YRR(3l.21).TSU3l.21),DAM,DAT 

♦ . YOSK 31,21 1.TS2! 31,21 ),WS131. 21 ).DEL( 31) 
common/ther5/wsmax(3l).tsnn(3l) 

DO 1000 I-l.NX 
ITER-0 

C TRANSFORMATION OF COORDINATE 

IF (KOU EQ.KST) THEN 
DO 100 J-l.NY 
YRR(I.J)=RT(J) 

GJUM. 

YTCJ) *Y0SUI.J)**.8/EXPIE/TS1(I.J)) 

100 TSRU)-TSl(I.J) 

CALL TRP2 (I.SUMY.DY. YT.NY.RIl.CJ) 

CALL TRPZ (I.SUMT.DY.TSR.NY.RIl.GJ) 

TTSUM(I)-OTS*RNK(I)*(TII)-TSR(1))+H2*SUMY*ROC1*SUMT+ROC2*TSR(NY) 
END IF 

B8-1./SQRTU. -DELHI) 

B82»8B*B8 
B83-(88*L J/tBB-L ) 

DO SCO J-2.NY 

IF (YOSUI.l) .EQ. 0.) THEN 

CHEC-OY 

ELSE 

CHEC-Y0SUU1/Y0SH 1. 1) 

END IF 

IF (CHEC .LE. .1JTHEN 

0ELT*(88*l. -IBB-l. )*BB3*11. -RT(J) ) )/(l. *B83**(1. -RT(J)) ) 

GO TO SS5 
ELSE 
DELT-L 
END IF 

SOO CONTINUE 
555 CONTINUE 

IF ((DELT ,LT. DEL(D). AND. (DELT .GE. DY)) DEL(I)-DELT 

IF ( DELT .LT. DY) DEUD-OY 

38-1./S0RTIL-0EU1)) 

382*88*88 
8B3H8BH. J/I88-L ) 

DO 200 J-l.NY 
YC*8B3**!l.-RTtJ)) 

YRR! I, J )-< BB+l. -< B8- 1. ) *YC)/( 1. * YC ) 

200 G J(J 1-2. *B8/(BB2-< 1. - YRRt I.J))**2. l/ALOGI BB3 ) 

TSS-TSl (I.D'Ll 

YSS»Y0(I)*.9S 

TSCT-O. 

TSLT-O. 

C 

C ITERATION ON ENERGY INTECRAL EQUATION 



100 




666 CONTINUE 

IF (ITER.EQ.30) CO TO 999 
YSGT-O. 

YSLT-O. 

ITER* (TER* l 
ITER1-0 

C ITERATION ON SPECIES EO. (SHOOT1NC METHOD) 



777 CONTINUE 
ITERl»lTERl‘l 

IF UTERI .GT. 30) CO TO 888 
BP-TSS-AP*YSS 
YOP»SH(I)*( YSS- YO( I ) )/G J( 1 ) 

CALL RUNCE (I.EYSR.AP.BP, YSS, Y0P,YP,NY.DY,BB.BB2.3B3.RII,G;,K0U) 

IF (A8S(YP).LT. .01) CO TO 888 

IF ((YP.LT. 0.). AND. UTERI .EQ.D) THEN 

YSS»Y0(I)*.99 

CO TO 777 

ENDIF 

IF (YP.CT.O.) THEN 

YSGT-YSS 

ELSE 

YSLT-YSS 

ENDIF 

YSS«(YSGT*YSLT)/2. 

CO TO 777 
888 CONTINUE 

C 

C END OF SPECIES EQUATION 



DO 300 J-l.NY 

300 TSR(J)»TSS-AP*(YSR(1)-YSR(J)) 

CALL TRP2 (I.SUMT.DY.TSR.NY.RI1.CJ) 
TNSUM«R0C1*SUMT»R0C2*TSR(NY) 

IF (ABS(TNSUM-TTSUM(U).LE LE-2) CO TO 999 
CORR-ABSITNSUM-TTSUM(I)) 

IF (CORR .LT. FIN) THEN 

FIN-CORR 

YFIN-YSS 

TFIN-TSS 

ENDIF 

IF (IITDLEQ.l) .AND. (TNSUM .LT. TTSUMII))) THEN 
TSS- TSS *2. 

CO TO 666 
ENDIF 

IF (TNSUM .CT. TTSUM(I)) THEN 

TSGT-TSS 

ELSE 

TSLT-TSS 

ENDIF 

IF ( (TSLT. EQ. 0. ). OR. (TSCT. EO. 0. ) ) THEN 

TNEW*TSS*.9S 

TNEW»<TSLT*TSGT)/2. 

ENDIF 

IF (ABS(TNEW-TSS) .LE t.E-3) CO TO 999 
IF (TNEW .CT. TSS) YSS-YSR(l) 

IF (TNEW .LT. TSS) YSS«YOSl(I.l)*.9 
TSS-TNEW 
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CO TO 666 
999 CONTINUE 
wmax*0. 

DO 400 J*l,NY 

WSUJ)- 8V*YSR( J ) •*. 8/EXP( E/T5R( J ) ) 

REF2- WS (U) 

WMAX- AMAX1 (REF2.WMAX) 

TSi (U)-TSRCJ) 

YOSi(U)-YSRU) 

400 YTCJ) *Y0S1(I*J) ## .8/EXP(E/T51(I,J)) 
do 410 j»l, nr 
410 ts2(i.j)*tsr(ny) 

CALL TRPZ (I.SUMY.DY, YT.NY.RII.GJ) 
TTSUM(I)»DTS*RNK(I)*{T(I)-TSR(i))*H2*SUMY«TTSUM(I) 

TS (I)-TSi(UJ 

YOS(I)-YOSUU) 

wsmax(i)»wmax 

C MAX SOLID TEMP IN RADIAL POSmON 

C 

call maxx (ny.tsr.km) 

TSRM(I)-RTUCM) 

tsmax(i)-tsrUcm) 

1000 CONTINUE 

C END OF ENERCY INTEGRAL EQUATION 

C MAX CAS and solid TEMP IN AXIAL POSITION 
C 

call maxx (nx,i,L) 
call maxx (nx. Umax. LI) 

C MAX SOLID and gas reaction rate IN AXIAL POSITION 



call maxx (nx,wsmax*m) 
call maxx (nx,wg,mi) 

RETURN 

END 

C INTEGRATION 

SUBROUTINE TRPZ (I.SUM.DY,FUN,N,RILGJ) 

DIMENSION FUN(2i),GJ(2I) 

COMMON/THER3/SH(3l),3C(31),YRR(3l,21).TSI(3L21 J.DAM.DAT 

♦ ,Y0Sl(3i,2I),TS2(3L21KWS(3L21),DEU3I) 

SUM«FUN(l)*RIl/2./GJU) 

DO 100 J-2.N-L 
R*<RI1+YRR(U))/GJ(J) 

100 SUM-SUM*FUN(J)*R 

SUM»iSUM^FUN(N) # (RH+YRR(I,N))/Z/GJ(N)) # DY 

RETURN 

END 

C RUNGE - KUTTA METHOD FOR SPECIES EQ. 

SUBROUTINE RUNGEI I, E. YSR. AP.BP. YSS, YOP, YP.NY, DR, BB, 9B2 BB3 RI1 GJ 
+ .KOU) 

DIMENSION YSR(2I).CJ(21) 

COMMON/THERi/ YH0I3I). Y0(3l). YHI3U. Y0S(3l). YH0V(3l). YH0SV(3t) 

♦ YICI3i),YU(3l),XJH(3l),RNK(3l),WC(31) t RO(31),TS(3n,Tt3l), 
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♦ XX(3U.DX(31),DY,DT,ROI.YHl, YOI,TI.CP.DS,XLS,XLDS.RT(3l) 
COUMON/TH£R3/SH(3l),BC(31).YRR(3l.21).TSH3l,2l).DAM,DAT 

♦ , YOSl(31.21).TS2(31.21), WS(3l.2I),DEL(31) 

F(Z)«Z 

G(A.B.Y.Z.TT)-(DAM*Y ,, .S/EXP(E/TT)-8*Z)/A/A 


YP-YOP 

YOSS-YSS 

YSRUJ-YSS 

DO 1000 J-l.NY-l 

RB-ORTLOATU-l) 

Cl-GJ(J) 

C2»Cl*U./(YRR(U)+Rn)-2.*a.-YRR(U!)/(BB2-(l.-YRR(U))**2.)) 

TT»AP‘YOSS«BP 

IF (YOSS.LT. 0. ) GO TO 999 

YI»F(YP)*DR 

Z1»G(C1,C2. YOSS.YP.TTJ'DR 
RBl»R8»0R/2. 

YC1»0B3**U.-RB1) 

RTl*(B8*l.-(B8-l. )*YC1)/(1. ♦YC1) 

a-2. •BB/(BB2-a.-RTl)"2. )/AL0G(BB3> 

C4»C3 , <1./IRT1*RI1)-2.*!I.-RT1)/(B82-(1.-RT1)**2.)) 

CHI*Y0SS*Yl/2- 

TT»AP‘GH1+BP 

IF (CH1.LT.0.) GO TO 999 

YZ»YP*Zl/2. 

Y2«F(YZ)*DR 

Z2»G(C3.C4,CH1,YZ,TT)*DR 

CH2«Y0SS*Y2/2. 

TT»AP*CH2*BP 

IF (CH2 .LT. 0.) GO TO 999 

YZ»YP*Z2/2. 

Y3»F(YZ!*DR 

Z3»G(C3.C4. CH2. YZ,TT!*DR 

CH3-Y0SS*Y3 

TT»AP*CH3*BP 

IF (CH3 .LT. 0.) GO TO 999 

YZ-YP-Z3 

Y4-F(YZ)*DR 

Z‘UG<a.C4.CH3,YZ,TT) , DR 

YOLD-YOSS 

YOSS- YOSS*( Y1 *2. *Y2*2. * Y3 * Y4 )/6. 

IF (YOSS .LT. 0.) GO TO 999 
IF (YOSS.GT. YOLD) THEN 
YP»l. 

GO TO 999 
END IF 

YP- YP*<Z1*2. *Z2+2. *Z3*Z4)/6. 

YSROH-YOSS 
1000 CONTINUE 
CO TO 777 
999 CONTINUE 

DO 888 K*J*l,NY 
388 YSROCMJ. 

777 CONTINUE 
RETURN 
END 
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TABLE 1 


Estimate of Transient Time Scales in a Monolithic Reactor 


Time 


Estimated 

Def inition Magni tude 


• * 


Gas residence time 

( i )( A . ) 

1-15 msec 

in reactor 

u A 

ref T 


Heat/mass transfer 
time between gas 

d* 2 


and solid surface 

a 

* 

0.5- 7.5 msec 

in channel 

Nu 7r4a (o,o) 

00 g 

d ^ 3 . 6 mm 

Mass diffusion time 

*2 

e h 


in porous catalytic 

c 

* 

0 . 03-0 . 3 msec 

layer 

D (o,o) 

e 

d £ 0 . 2 mm 

s 

Heat conduction time 

h* 2 


in porous catalytic 

c 

# 

0.1- 4 msec 

layer 

a 

c 


Heat conduction time 

h‘ 2 


in substrate layer 

s 

* 

0.4- 21 msec 


a 

s 

Same order as 

Homogeneous chemical 


gas residence 

reaction time 


time 


(continued on next page) 
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Time 


Heterogeneous 
reaction time 


TABLE 1 (continued) 


Definition 


Estimated 

Magnitude 


Same order as 

chemical solid heat-up 

time 


2 



* * * * 

p C r dr 


r 

d 


♦ 

Nu A (o,o) 
oo q 


Solid heat-up time 


0.5- 20 sec 
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TABLE 2 


List of 

Time Scales 

and Nondimensional 



Parameters 

in Case 1 


Name 

Symbol 

Value 

Unit 

Mass diffusion time 
in catalytic layer 

* 

T 

m 

4.3xl0" 5 

sec 

Heat conduction time 
in catalytic layer 

* 

T 

t 

4.3xl0 -5 

sec 

Heat conduction time 
in substrate layer 

* 

T 

s 

N> 

X 

O 

i 

•t* 

sec 

Gas residence time 

* 

T 

r 

3. xl0~ 3 

sec 

Solid heat-up time 

* 

T 

HT 

2 . x 10~ 2 

sec 

Damkohler number 
of gas 

D 

am 

0.04 


Damkohler number 
of solid 

D 

at 

0.04 


Sherwood number 

Sh 

2.36 - 0.27 


Biot number 

Be 

o 

v£> 

-sj 

* 

O 
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TABLE 3 


Parametric Studies 


Parameter 

Case 1 

Case 2 

Case 3 

Case 4 

Case 5 

0*(kg-Pt/m 3 ) 

1.8 

0. 18 

18. 

1.8 

1.8 

C* (cal/kg/°K) 

1850. 

1850. 

1850. 

1310. 

1850. 

P 

T* (°K) 

1372. 

1372. 

1372. 

2200. 

1372 

ab 

* 

r (sec) 

0.02 

0.02 

0.02 

0.02 

2. 

HT 

<P 

8. 

8. 

8. 

4. 

8. 

Dam 

0.04 

0.004 

0.4 

0.04 

0.04 

Dat 

0.04 

0.004 

0.4 

0.08 

0.04 

q 

c 

5. 

5. 

5. 

7. 

5. 

Y (o,o) 

0.5 

0.5 

0.5 

0.66 

0.5 

g 

f 

D 

0.08 

0.008 

0.8 

0.06 

0.08 

am 

/ 

D 

0.2 

0.02 

2. 

0.56 

0.2 

a t 

Max. Sh 

2.36 

2.36 

2.36 

2.8 

2.36 

Min. Sh 

0.27 

0.27 

0.27 

0.27 

0.27 

Max. Be 

0.97 

0.97 

0.97 

1 . 15 

0.97 

Min. Be 

0. 11 

0. 11 

0. 11 

0. 11 

0. 11 
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TABLE 4 


List of Property Values in Case 1 


Name 

Symbol 

Value 

Unit 

Diameter of cell 

* 

d 

10- 3 

m 

Length of cell 

* 

L 

5x10" 2 

m 

Porosity of catalytic layer 

€ 

0.5 


Thickness of catalytic layer 

* 

h 

c 

3. 3xl0~ 3 

m 

Thickness of substrate layer 

* 

h 

s 

ON 

-s) 

X 

o 

1 

u 

m 

Open area ratio 

« * 

A /A 

s 

2. 


Density of catalytic layer 

* 

p 

c 

13. 

kg/m 3 

Density of substrate layer 

* 

p 

5 

183. 

kg/m 3 

Specific heat of catalytic layer 

* 

c 

c 

180. 

cal/kg/°K 

Specific heat of substrate layer 

* 

c 

s 

180. 

cal/kg/°K 

Heat conductivity of 
catalytic layer 

* 

A 

c 

0.069 

cal/m/s/°K 

Heat conductivity of 
substrate layer 

* 

A 

s 

0.69 

cal/m/s/°K 

Effective diffusivity 
of oxygen in catalytic 
layer 

D (o,o) 

e 

1 . 66x 10” 5 

m 2 /sec 
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TABLE 5 


List of Operating Parameters in Case 1 


Name 

Symbol 

Value 

Unit 

Catalyst loading 

* 

P 

1.8 

kg/m 3 

Inlet pressure 

P*(0,0) 

9. 

atm 

Inlet Temperature 

T*(o,o) 

400. 

°K 

Inlet velocity 

* 

U (0,0) 

10. 

m/sec 

Inlet mass ratio 


1 . 


Hydrogen/Oxygen 
equivalence ratio 

4 > 

8. 


Adiabatic flame 
temperature 

* 

T 

ab 

1376 

°K 
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TABLE A. 1 


Parameters 

for Two Sets 

of Gas-Phase 

Reaction Rates 

Parameter 

Symbol 

Value 

Unit 

Rate Constant 

• 

K 

L 

CO 

o 

T — ( 

X 

atm 1 ‘ 5 (kg/m 3 )/: 


» 

K 

H 

vO 

X 

* — * 
o 

o\ 

1/atm/sec 

Exponent of P 

a 

L 

-1.5 



a 

H 

1 


* 

Exponent of p 

0 

b 

L 

0 . 1 



b 

H 

1 . 


Activation Energy 

• 

E 

L 

7 . 9x10* 

cal/mole 


* 

E 

H 

1.9xl0 4 

cal/mole 


Note : 

Subscript L and H denote the parameter used for the gas-phase 

* * » * 

reaction rate at T ^ T and T ^ T respectively. The critical 

cr cr 

* Q 

temperature in this study is T = 1200 K. The general expression 

cr 

for the gas-phase reaction rate is : 

W* = K* p* a P* b exp ( -E*/R*T* ) [kg-0 /m 3 /sec] 

o 2 




110 




gas channel 


sol id 
substrate 

catalytic layer 


(a) Physical model for a tubular reactor 



(b) Illustration of catalytic layer structure 


gas channel 


] 

r 

d i 


1 

catalytic layer i 


substrate layer 




t 






center of 
gas channel 


(c) Defination of dimensions 


Figure 1 : Monolithic reactor model description 














Figure 4 : 


Solid contours for Case 1 at t 








m t n 
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c 

0 


c 

o 


c 



Catalyst loading (0* kg-pt/* 3 ) 


Figure 18 : Ignition boundary map 

( u « 10 m/sec) 

CD 


a : Gas-chcse ignition 
C : Ca*a:y:ic igniiicr: 
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Figure 23 : Solig in-depjh temperature profiles for Case 3 

at t = 0.2 t 

HT 

Solid line (x=0), « (x=0.005), o(x=0.01), A (x=0.5) 







r 


Figure 25 : Solitj in-d|pth temperature profiles for Case 3 

at t “It 

ht 

Solid line (x=0), • (x=0.005), o (x=0.01), (x=0.5) 



Figure 26 : Soli^ in-d|pth species profiles for Case 3 

at t = 1 t 

HT 

Solid line (x=0), • (x=0.005), « (x=0.01) 




r 


Figure 27 : Soli$ in-dfpth temperature profiles for Case 3 

at 1 = 3 T ht a 

Solid line (x=0), • (x=0.005), •(x-0.01), Q (x=0.5) 







r 


Figure 29 : SoliQ in-dep^h temperature profiles for Case 4 

at t = 0. 2 x 

H T 

Solid line (x=0), • (x=0.06), o (x=0. 12) , A (x=0.S) 




O.C 0.2 0.4 0.6 0.8 

r 


Figure 30 : Soli^ in-depjh species profiles for Case 4 

at t ■ 0.2 x 

HT 

Solid line (x=0), • (x=0.06), o (x=0. 12) 






0.0 1.0 2,0 3.0 


r 


Figure 31 : Solid # in-dejth temperature profiles f or Case 4 

at t ■ 1 r 

HT 

Solid line (x=0), • (x=0.06), o (x=0.12), (x=0.5) 
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SolicJ in-depth temperature profiles for Case 4 
at 1 = 2 T ht 

Solid line (x=0), • (x*0.06), o (x=0.12), Q (x=0.5) 


Figure 33 : 
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Figure 34 : Soli^ in-djpth species profiles for Case 4 

at t = 2 t 

HT 

Solid line (x=0), • (x=0.06), o (x=0.5) 



X 


Figure 35 : Comparisons of tj/o models for Case 1 on T 

at t - 1. 2, 3 t n 5 

Solid line : fir! transient mocel 
Symbol O : energy integral mode! 



Figure 36 : Comparisons of tyo models for Case 1 on y 

at t » 1, 2, 3 r g 

HT 

Solid line : full transient mocel 
Symbol O : energy intec r cl mode 




0.0 u.z U . u u.o u. o 

X 


Figure 37 : Comparisons of tyo models for Case 

at t = 1 , 2, 3 


.u 


on W 

<3 


Solid line : full transient model 
Symbol O : energy integral mode 
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Figure 39 : Comparisons of two models # for C^se 1 on solid 

in-depth temperature at t * 1 t ht 

Solid line : full transient model 
Symbol o : energy integral model 




Figure 40 : Comparisons of two models f^r Case 1 on solid 

ln-depth species at t * 1 t ht 

Solid line : full transient model 
Symbol o : energy integral model 
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Figure 42 : Comparisons of two models for Case 1 on solid 

in-depth temperature at t = 2 x* 

HT 

Solid line : full transient model 
Symbol o : energy integral model 





Figure 43 : Comparisons of two models f or # Case 1 on solid 

in-depth species at t * 2 t ht 

Solid line : full transient model 
Symbol o : energy integral model 
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Figure 45 : Comparisons of two models for Case 1 on solid 

in-depth temperature at t = 3 t ht 

Solid line : full transient model 
Symbol o : energy integral model 
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r 


Figure 46 : Comparisons of two models fop Case 1 on solid 

in-depth species at t = 3 t* 

ht 

Solid line : full transient model 
Symbol o : energy integral model 



0.0 


0.2 


0 .^ 


0.6 0.8 

X 


Figure 47 : Comparisons of two # models for Case 

at t = 0.2, 1, 3 t ht 

Sol id line : full transient model 
Symbol o : energy integral model 


on T 

t 




Figure 48 : 

Comparisons 

of two Bodels for Case 3 on y 

• < 


at t 3 0.2, 

l ’ 3 t ht 


Solid line : 

: full transient model 


Symbol o 

: energy integral model 




0.0 0.2 0 A 0.6 0.8 
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Figure 49 : Comparisons of two^models for Case 3 on W 

at t = 0.2, 1, 3 t 

HT 

Solid line : full transient model 
Symbol o : energy Integral model 









3.C 
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Figure 51 : Comparisons of two models^for Cas£ 3 on solid 

in-depth temperature at t = 0.2 t ht 

Solid line : full transient model 
Symbol o : energy integral model 











163 



r 


Figure 54 : Comparisons of two models for Ca sf 3 on solid 

in-depth temperature at t = 1 t ht 

Solid line : full transient model 
Symbol o : energy integral model 




Figure 55 : Comparisons of two models fop Case 3 on solid 

in-depth species at t = 1 x 

H T 

Solid line : full transient model 
Symbol o : energy integral model 
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0.0 0.2 0.4 0.6 0.3 

r 


Figure 58 : Comparisons of two models for Case 3 on solid 

in-depth species at t *3 t ht 

Solid line : full transient model 
Symbol o : energy integral model 












(3) 1/T (K) 


(3) 1/T (K) 


Figure A.l : Comparisons of ignition delay time between model 

predicted by Eqs. (A. 10), (A. 11) and experimental 
data [14], [IS]. 

Solid line : Model prediction 

(A) 47. H2, 27.02, o : p » 1. • : p « 2 atm 

(B) 8 7. H2, 27.02, o:p-l. • : p « 5 atm 

(C) 127. H 2 , 27.02, o:p«l, • : p ■ 2 atm 

(D) 677. H2, 337.02, • : p - 1. • : p « 2 atm 
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